{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exploratury Data Analysis - Part II\n",
    "By Minaduki 2019.    \n",
    "项目的Git存储库：https://git.nju.edu.cn/Minaduki/datascience     \n",
    "本作业数据来源自：https://github.com/cmawer/pycon-2017-eda-tutorial       \n",
    "This is a EDA using the Food and Agriculture Organization(FAO) of the United Nation's AQUASTAT dataset.    \n",
    "The data used can be found [here](http://www.fao.org/nr/water/aquastat/data/query/index.html).      "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import stats\n",
    "from scipy.stats import t\n",
    "import math\n",
    "import statistics\n",
    "import random"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 准备工作\n",
    "导入AQUASTAT数据，并进行如下几步操作：\n",
    "1. 将表中的几个统计的时间节点保存到一个列表中，方便之后调用和随机操作。\n",
    "2. 使用上次作业中编写的数据切分函数，将数据表中的人均GDP字段的值切分出来，存储到新的数据表中。\n",
    "3. 为了方便起见，将数据表中含有空值的条目全部丢弃。\n",
    "4. 显示新的数据表的信息和一些数据，检查导入是否成功。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_data = pd.read_csv('aquastat.csv.gzip', compression='gzip')\n",
    "time_periods = my_data.time_period.unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def slice_by_variable(dataframe, variable):\n",
    "    dataframe = dataframe[dataframe.variable==variable]\n",
    "    dataframe = dataframe.pivot(index='country', columns='time_period', values='value')\n",
    "    return dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def slice_by_country(dataframe, country):\n",
    "    dataframe = dataframe[dataframe.country==country] \n",
    "    dataframe = dataframe.pivot(index='variable', columns='time_period', values='value')\n",
    "    dataframe.index.name = country\n",
    "    return dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def slice_by_time(dataframe, time_period):\n",
    "    dataframe = dataframe[dataframe.time_period==time_period] \n",
    "    dataframe = dataframe.pivot(index='country', columns='variable', values='value')\n",
    "    dataframe.columns.name = time_period\n",
    "    return dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'pandas.core.frame.DataFrame'>\n",
      "Index: 108 entries, Algeria to Zimbabwe\n",
      "Data columns (total 12 columns):\n",
      "1958-1962    108 non-null float64\n",
      "1963-1967    108 non-null float64\n",
      "1968-1972    108 non-null float64\n",
      "1973-1977    108 non-null float64\n",
      "1978-1982    108 non-null float64\n",
      "1983-1987    108 non-null float64\n",
      "1988-1992    108 non-null float64\n",
      "1993-1997    108 non-null float64\n",
      "1998-2002    108 non-null float64\n",
      "2003-2007    108 non-null float64\n",
      "2008-2012    108 non-null float64\n",
      "2013-2017    108 non-null float64\n",
      "dtypes: float64(12)\n",
      "memory usage: 10.5+ KB\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>time_period</th>\n",
       "      <th>1958-1962</th>\n",
       "      <th>1963-1967</th>\n",
       "      <th>1968-1972</th>\n",
       "      <th>1973-1977</th>\n",
       "      <th>1978-1982</th>\n",
       "      <th>1983-1987</th>\n",
       "      <th>1988-1992</th>\n",
       "      <th>1993-1997</th>\n",
       "      <th>1998-2002</th>\n",
       "      <th>2003-2007</th>\n",
       "      <th>2008-2012</th>\n",
       "      <th>2013-2017</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>country</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>Algeria</td>\n",
       "      <td>171.20</td>\n",
       "      <td>252.4</td>\n",
       "      <td>439.70</td>\n",
       "      <td>1186.0</td>\n",
       "      <td>2197.0</td>\n",
       "      <td>2790.0</td>\n",
       "      <td>1766.0</td>\n",
       "      <td>1612.0</td>\n",
       "      <td>1774.0</td>\n",
       "      <td>3940.0</td>\n",
       "      <td>5582.0</td>\n",
       "      <td>4210.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>Argentina</td>\n",
       "      <td>1149.00</td>\n",
       "      <td>1058.0</td>\n",
       "      <td>1401.00</td>\n",
       "      <td>2112.0</td>\n",
       "      <td>2908.0</td>\n",
       "      <td>3543.0</td>\n",
       "      <td>6804.0</td>\n",
       "      <td>8177.0</td>\n",
       "      <td>2579.0</td>\n",
       "      <td>8231.0</td>\n",
       "      <td>14348.0</td>\n",
       "      <td>12622.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>Australia</td>\n",
       "      <td>1860.00</td>\n",
       "      <td>2537.0</td>\n",
       "      <td>3886.00</td>\n",
       "      <td>7740.0</td>\n",
       "      <td>12847.0</td>\n",
       "      <td>11584.0</td>\n",
       "      <td>18531.0</td>\n",
       "      <td>23551.0</td>\n",
       "      <td>20191.0</td>\n",
       "      <td>40666.0</td>\n",
       "      <td>67217.0</td>\n",
       "      <td>55906.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>Austria</td>\n",
       "      <td>1085.00</td>\n",
       "      <td>1567.0</td>\n",
       "      <td>2906.00</td>\n",
       "      <td>6743.0</td>\n",
       "      <td>9369.0</td>\n",
       "      <td>16264.0</td>\n",
       "      <td>24958.0</td>\n",
       "      <td>26450.0</td>\n",
       "      <td>26248.0</td>\n",
       "      <td>46500.0</td>\n",
       "      <td>48137.0</td>\n",
       "      <td>43768.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>Bahamas</td>\n",
       "      <td>1753.00</td>\n",
       "      <td>2557.0</td>\n",
       "      <td>3323.00</td>\n",
       "      <td>3617.0</td>\n",
       "      <td>7164.0</td>\n",
       "      <td>11150.0</td>\n",
       "      <td>11684.0</td>\n",
       "      <td>17286.0</td>\n",
       "      <td>22503.0</td>\n",
       "      <td>24303.0</td>\n",
       "      <td>22112.0</td>\n",
       "      <td>22898.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>Bangladesh</td>\n",
       "      <td>99.58</td>\n",
       "      <td>121.7</td>\n",
       "      <td>93.05</td>\n",
       "      <td>128.9</td>\n",
       "      <td>215.7</td>\n",
       "      <td>247.6</td>\n",
       "      <td>285.7</td>\n",
       "      <td>390.4</td>\n",
       "      <td>401.7</td>\n",
       "      <td>543.1</td>\n",
       "      <td>856.6</td>\n",
       "      <td>1211.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>Barbados</td>\n",
       "      <td>416.80</td>\n",
       "      <td>579.5</td>\n",
       "      <td>966.80</td>\n",
       "      <td>1990.0</td>\n",
       "      <td>4586.0</td>\n",
       "      <td>6616.0</td>\n",
       "      <td>7464.0</td>\n",
       "      <td>9551.0</td>\n",
       "      <td>11674.0</td>\n",
       "      <td>16459.0</td>\n",
       "      <td>15316.0</td>\n",
       "      <td>15662.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>Belgium</td>\n",
       "      <td>1443.00</td>\n",
       "      <td>2110.0</td>\n",
       "      <td>3876.00</td>\n",
       "      <td>8540.0</td>\n",
       "      <td>9428.0</td>\n",
       "      <td>15201.0</td>\n",
       "      <td>23562.0</td>\n",
       "      <td>25017.0</td>\n",
       "      <td>24988.0</td>\n",
       "      <td>44092.0</td>\n",
       "      <td>44946.0</td>\n",
       "      <td>40181.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>Belize</td>\n",
       "      <td>327.10</td>\n",
       "      <td>420.4</td>\n",
       "      <td>519.40</td>\n",
       "      <td>858.8</td>\n",
       "      <td>1183.0</td>\n",
       "      <td>1587.0</td>\n",
       "      <td>2667.0</td>\n",
       "      <td>2953.0</td>\n",
       "      <td>3557.0</td>\n",
       "      <td>4325.0</td>\n",
       "      <td>4674.0</td>\n",
       "      <td>4907.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>Benin</td>\n",
       "      <td>94.46</td>\n",
       "      <td>112.0</td>\n",
       "      <td>135.00</td>\n",
       "      <td>218.5</td>\n",
       "      <td>322.5</td>\n",
       "      <td>343.7</td>\n",
       "      <td>315.2</td>\n",
       "      <td>356.6</td>\n",
       "      <td>411.9</td>\n",
       "      <td>685.5</td>\n",
       "      <td>807.7</td>\n",
       "      <td>779.1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "time_period  1958-1962  1963-1967  1968-1972  1973-1977  1978-1982  1983-1987  \\\n",
       "country                                                                         \n",
       "Algeria         171.20      252.4     439.70     1186.0     2197.0     2790.0   \n",
       "Argentina      1149.00     1058.0    1401.00     2112.0     2908.0     3543.0   \n",
       "Australia      1860.00     2537.0    3886.00     7740.0    12847.0    11584.0   \n",
       "Austria        1085.00     1567.0    2906.00     6743.0     9369.0    16264.0   \n",
       "Bahamas        1753.00     2557.0    3323.00     3617.0     7164.0    11150.0   \n",
       "Bangladesh       99.58      121.7      93.05      128.9      215.7      247.6   \n",
       "Barbados        416.80      579.5     966.80     1990.0     4586.0     6616.0   \n",
       "Belgium        1443.00     2110.0    3876.00     8540.0     9428.0    15201.0   \n",
       "Belize          327.10      420.4     519.40      858.8     1183.0     1587.0   \n",
       "Benin            94.46      112.0     135.00      218.5      322.5      343.7   \n",
       "\n",
       "time_period  1988-1992  1993-1997  1998-2002  2003-2007  2008-2012  2013-2017  \n",
       "country                                                                        \n",
       "Algeria         1766.0     1612.0     1774.0     3940.0     5582.0     4210.0  \n",
       "Argentina       6804.0     8177.0     2579.0     8231.0    14348.0    12622.0  \n",
       "Australia      18531.0    23551.0    20191.0    40666.0    67217.0    55906.0  \n",
       "Austria        24958.0    26450.0    26248.0    46500.0    48137.0    43768.0  \n",
       "Bahamas        11684.0    17286.0    22503.0    24303.0    22112.0    22898.0  \n",
       "Bangladesh       285.7      390.4      401.7      543.1      856.6     1211.0  \n",
       "Barbados        7464.0     9551.0    11674.0    16459.0    15316.0    15662.0  \n",
       "Belgium        23562.0    25017.0    24988.0    44092.0    44946.0    40181.0  \n",
       "Belize          2667.0     2953.0     3557.0     4325.0     4674.0     4907.0  \n",
       "Benin            315.2      356.6      411.9      685.5      807.7      779.1  "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_gdp = slice_by_variable(my_data, 'gdp_per_capita')\n",
    "data_gdp = data_gdp.dropna(axis=0)\n",
    "data_gdp.info()\n",
    "data_gdp.head(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 抽取数据\n",
    "随机抽取一个年份的数据存入变量data中，作为之后使用的数据，打印出数据的概况和位置性、离散性测度信息如下所示："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[1;31m统计区间： 1998-2002 \u001b[0m\n",
      "country\n",
      "Algeria                                1774.0\n",
      "Argentina                              2579.0\n",
      "Australia                             20191.0\n",
      "Austria                               26248.0\n",
      "Bahamas                               22503.0\n",
      "                                       ...   \n",
      "United States of America              38132.0\n",
      "Uruguay                                4088.0\n",
      "Venezuela (Bolivarian Republic of)     3657.0\n",
      "Zambia                                  376.5\n",
      "Zimbabwe                                499.7\n",
      "Name: 1998-2002, Length: 108, dtype: float64\n",
      "\n",
      "\u001b[1;31m位置性测度统计结果：\u001b[0m\n",
      "均值：\t\t 8421.026851851851\n",
      "中位数：\t 2532.0\n",
      "第25个百分位数： 526.3000000000001\n",
      "\n",
      "\u001b[1;31m离散性测度统计结果：\u001b[0m\n",
      "变化范围：\t [ 115.3 \t 52532.0 ]\n",
      "极差：\t\t 52416.7\n",
      "方差：\t\t 138308703.7840386\n",
      "标准差：\t 11760.472090185776\n",
      "变异系数：\t 1.396560336059201\n"
     ]
    }
   ],
   "source": [
    "temp_year = random.randint(0,11)\n",
    "print('\\033[1;31m统计区间：',time_periods[temp_year],'\\033[0m')\n",
    "data = data_gdp[time_periods[temp_year]]\n",
    "print(data)\n",
    "print('\\n\\033[1;31m位置性测度统计结果：\\033[0m')\n",
    "print('均值：\\t\\t',data.mean())\n",
    "print('中位数：\\t',data.median())\n",
    "print('第25个百分位数：',data.quantile(q=0.25))\n",
    "print('\\n\\033[1;31m离散性测度统计结果：\\033[0m')\n",
    "print('变化范围：\\t [',data.min(),'\\t',data.max(),']')\n",
    "print('极差：\\t\\t',data.max()-data.min())\n",
    "print('方差：\\t\\t',data.var())\n",
    "print('标准差：\\t',data.std())\n",
    "print('变异系数：\\t',data.std()/data.mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 随机数生成\n",
    "接下来做一些均值估计的准备工作，首先指定一个随机数种子，然后先生成一个容量为80的小样本："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(4396)\n",
    "sample_data=np.random.choice(a=data,size=50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "打印出生成的样本的均值看看："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "样本的均值为:  7840.018\n"
     ]
    }
   ],
   "source": [
    "print('样本的均值为: ',sample_data.mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "再看看多个样本的均值分布，这里与课上一样，也取500次样本："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "样本均值的均值为:  8375.456167999999\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAFSpJREFUeJzt3X+M5HV9x/HnW0A9Wb0Dgc150p5GaqRcRG9DQROzqxZRG3802kCNQtWc9Ve0vbRiTVOtbYL1V9LWVE9RiaGsVqUQ/FWKbKlNq91T4O48EdCrcpx3UvBkqUEP3/1jvouzy87O7Mz3uzP7uecjmezM5/ud2dfNzr7uu9/vZ74TmYkkae172LADSJLqYaFLUiEsdEkqhIUuSYWw0CWpEBa6JBXCQpekQljoOqpFxIkRcWVE3BcR/xMRvz/sTFK/jh12AGnIPgT8HBgHzgS+EBE3Zeae4caSVi58p6iOVhFxPHAPcEZmfrca+xSwPzMvHmo4qQ/uctHR7DeAB+bLvHIT8JtDyiMNxELX0WwMOLxo7DDw6CFkkQZmoetoNgc8ZtHYY4B7h5BFGpiFrqPZd4FjI+K0trGnAh4Q1ZrkQVEd1SJiGkjgtbRmuXwReIazXLQWuYWuo90bgHXAIeAK4PWWudYqt9AlqRBuoUtSISx0SSqEhS5JhbDQJakQq3pyrpNOOik3b968YOy+++7j+OOPX80YfTNrM8zaDLPWb1g5d+7ceVdmntx1xcxctcvWrVtzseuvv/4hY6PKrM0wazPMWr9h5QRms4eOdZeLJBXCQpekQljoklQIC12SCmGhS1IhLHRJKkTXQo+IR0bENyLipojYExHvqsY/GRHfj4gbq8uZzceVJHXSyxuL7geenZlzEXEc8LWI+FK17E8y87PNxZMk9aproVeT2ueqm8dVF8+5K0kjpqfzoUfEMcBO4EnAhzLzbRHxSeAcWlvw1wEXZ+b9S9x3G7ANYHx8fOv09PSC5XNzc4yNjQ34z1gdZm3GSrLu2r/4M51/Zcum9XVF6qjU53XY1krWYeWcmpramZkT3dZb0QdcRMQG4ErgzcD/Aj8CHg7sAG7PzL9c7v4TExM5Ozu7YGxmZobJycmeMwyTWZuxkqybL/5Cx2X7LnlhTYk6K/V5Hba1knVYOSOip0Jf0SyXzPwJMAOcl5kHqtMM3A98Ajirr6SSpFr0Msvl5GrLnIhYBzwX+E5EbKzGAngJsLvJoJKk5fUyy2UjcFm1H/1hwGcy85qI+GpEnAwEcCPwhw3mlCR10cssl5uBpy0x/uxGEkmS+uI7RSWpEBa6JBXCQpekQljoklQIC12SCmGhS1IhepmHLq2a5d7aL2l5bqFLUiEsdEkqhIUuSYWw0CWpEBa6JBXCQpekQljoklQIC12SCmGhS1IhLHRJKoSFLkmFsNAlqRAWuiQVomuhR8QjI+IbEXFTROyJiHdV40+IiK9HxK0R8emIeHjzcSVJnfSyhX4/8OzMfCpwJnBeRJwNvAf4YGaeBtwDvKa5mJKkbroWerbMVTePqy4JPBv4bDV+GfCSRhJKknoSmdl9pYhjgJ3Ak4APAe8F/iszn1QtPxX4UmaescR9twHbAMbHx7dOT08vWD43N8fY2NiA/4zVYdZmtGfdtf9w34+zZdP6uiJ1tFaf11G3VrIOK+fU1NTOzJzotl5Pn1iUmQ8AZ0bEBuBK4ClLrdbhvjuAHQATExM5OTm5YPnMzAyLx0aVWZvRnvWiAT6xaN8rJusJtIy1+ryOurWSddRzrmiWS2b+BJgBzgY2RMT8fwiPB+6sN5okaSV6meVycrVlTkSsA54L7AWuB15WrXYhcFVTISVJ3fWyy2UjcFm1H/1hwGcy85qI+DYwHRF/BXwLuLTBnJKkLroWembeDDxtifHvAWc1EUqStHI9HRSVSrB5mQOu+y554SomkZrhW/8lqRAWuiQVwkKXpEJY6JJUCAtdkgrhLBepBs6g0ShwC12SCmGhS1IhLHRJKoSFLkmFsNAlqRDOclHtnPEhDYdb6JJUCAtdkgphoUtSISx0SSqEhS5JhbDQJakQFrokFcJCl6RCdC30iDg1Iq6PiL0RsSci3lKNvzMi9kfEjdXlBc3HlSR10ss7RY8A2zPzmxHxaGBnRFxbLftgZr6vuXiSpF51LfTMPAAcqK7fGxF7gU1NB5MkrUxkZu8rR2wGbgDOAP4YuAj4KTBLayv+niXusw3YBjA+Pr51enp6wfK5uTnGxsb6Cr/azNqbXfsPd1y2ZdP6h4y1Z13uvt0s9diD5FpKp+e1jseum6/X+g0r59TU1M7MnOi2Xs+FHhFjwL8Bf52Zn4+IceAuIIF3Axsz89XLPcbExETOzs4uGJuZmWFycrKnDMNm1t6s9ORc7VmXu2833U78VcdJwzo9r6N4QjJfr/UbVs6I6KnQe5rlEhHHAZ8DLs/MzwNk5sHMfCAzfwl8FDhrkMCSpMH0MsslgEuBvZn5gbbxjW2rvRTYXX88SVKvepnl8kzglcCuiLixGvsz4IKIOJPWLpd9wOsaSShJ6kkvs1y+BsQSi75YfxxJUr/8xCKt2CAHLps0SK5u9/WTlrQW+NZ/SSqEhS5JhbDQJakQFrokFcJCl6RCWOiSVAgLXZIKYaFLUiEsdEkqhIUuSYWw0CWpEJ7LRerB/Lletm85wkUjei4byS10SSqEhS5JhbDQJakQFrokFcJCl6RCOMvlKLXcJ/Q0+ek8S31fZ45I9XALXZIK0bXQI+LUiLg+IvZGxJ6IeEs1fmJEXBsRt1ZfT2g+riSpk1620I8A2zPzKcDZwBsj4nTgYuC6zDwNuK66LUkakq6FnpkHMvOb1fV7gb3AJuDFwGXVapcBL2kqpCSpu8jM3leO2AzcAJwB/CAzN7QtuyczH7LbJSK2AdsAxsfHt05PTy9YPjc3x9jYWD/ZV91aynro7sMc/Fl/992yaf2yy3ftP9zfA3cwvo6+s662frJ2ez6bspZer2sl67ByTk1N7czMiW7r9TzLJSLGgM8Bb83Mn0ZET/fLzB3ADoCJiYmcnJxcsHxmZobFY6NqLWX9u8uv4v27+pvEtO8Vk8sur3tGyvYtR/rOutr6ydrt+WzKWnq9rpWso56zp1kuEXEcrTK/PDM/Xw0fjIiN1fKNwKFmIkqSetHLLJcALgX2ZuYH2hZdDVxYXb8QuKr+eJKkXvXyt+MzgVcCuyLixmrsz4BLgM9ExGuAHwAvbyaiJKkXXQs9M78GdNph/px640iS+rU2jkRpVS13WgCtXLfns8lTLejo4lv/JakQFrokFcJCl6RCWOiSVAgLXZIK4SwXaciG9WEjKo9b6JJUCAtdkgphoUtSISx0SSqEhS5JhXCWyxB5jg9JdXILXZIKYaFLUiEsdEkqhIUuSYWw0CWpEM5ykUaYM6G0Em6hS1IhLHRJKkTXQo+Ij0fEoYjY3Tb2zojYHxE3VpcXNBtTktRNL1vonwTOW2L8g5l5ZnX5Yr2xJEkr1bXQM/MG4O5VyCJJGkBkZveVIjYD12TmGdXtdwIXAT8FZoHtmXlPh/tuA7YBjI+Pb52enl6wfG5ujrGxsX7zr6q6s+7af3jZ5Vs2re/7sQ/dfZiDP+v77qtqfB1m7dNyr5Gj+XerKcPKOTU1tTMzJ7qt12+hjwN3AQm8G9iYma/u9jgTExM5Ozu7YGxmZobJycmuGUZB3VmbnJL2d5dfxft3rY1Zqdu3HDFrn5Z7jRzNv1tNGVbOiOip0Pua5ZKZBzPzgcz8JfBR4Kx+HkeSVJ++Cj0iNrbdfCmwu9O6kqTV0fVvx4i4ApgEToqIO4C/ACYj4kxau1z2Aa9rMKMkqQddCz0zL1hi+NIGskiSBjA6R3e0YssdVN2+ZRWDSBoJvvVfkgphoUtSISx0SSqEhS5JhbDQJakQFrokFcJCl6RCWOiSVAgLXZIKYaFLUiEsdEkqhOdykdaw5c/nc4TJ1YuiEeAWuiQVwkKXpEJY6JJUCAtdkgphoUtSISx0SSqEhS5Jheha6BHx8Yg4FBG728ZOjIhrI+LW6usJzcaUJHXTyxb6J4HzFo1dDFyXmacB11W3JUlD1LXQM/MG4O5Fwy8GLquuXwa8pOZckqQV6ncf+nhmHgCovp5SXyRJUj8iM7uvFLEZuCYzz6hu/yQzN7Qtvyczl9yPHhHbgG0A4+PjW6enpxcsn5ubY2xsrN/8q6rurLv2H67tsRYbXwcHf9bYw9fKrM0YXwennLh+2DF6slZ6YFg5p6amdmbmRLf1+j0518GI2JiZByJiI3Co04qZuQPYATAxMZGTk5MLls/MzLB4bFTVnfWiZU6sNKjtW47w/l1r49xrZm3G9i1H+L2j9HerKaOes99dLlcDF1bXLwSuqieOJKlfvUxbvAL4T+DJEXFHRLwGuAT47Yi4Ffjt6rYkaYi6/u2YmRd0WPScmrNIkgbgO0UlqRBr4+iOpKIs/qSl7VuOPDhJYN8lLxxGpCK4hS5JhbDQJakQFrokFcJCl6RCeFBUKtjig4/tPPhYHrfQJakQFrokFcJCl6RCWOiSVAgLXZIK4SwX6Si13AwYcBbMWuQWuiQVwkKXpEJY6JJUCAtdkgphoUtSIZzl0oP52QDtJ+Gf120mQLeZBNJa5AyZ0eQWuiQVwkKXpEIMtMslIvYB9wIPAEcyc6KOUJKklatjH/pUZt5Vw+NIkgbgLhdJKkRkZv93jvg+cA+QwEcyc8cS62wDtgGMj49vnZ6eXrB8bm6OsbGxvjOshl37DwMwvg4O/mzhsi2b1vd039W2VNZRZdZmDJp1udd2t9f1Sn8v2rN2u+8wDauvpqamdvayS3vQQn9cZt4ZEacA1wJvzswbOq0/MTGRs7OzC8ZmZmaYnJzsO8NqaJ+2+P5dC/dSjeq0xaWyjiqzNmPQrMu9tgedtrj4/u1ZR3nK47D6KiJ6KvSBdrlk5p3V10PAlcBZgzyeJKl/fRd6RBwfEY+evw6cC+yuK5gkaWUG+dtxHLgyIuYf5x8z88u1pJIkrVjfhZ6Z3wOeWmMWSdIA1sbRnR4sd5BmlA+ySKPK8xCtPc5Dl6RCWOiSVAgLXZIKYaFLUiEsdEkqxJqZ5dLkEXeP5kv18ndqONxCl6RCWOiSVAgLXZIKYaFLUiEsdEkqxJqZ5TKqPJov1WvQD884mrmFLkmFsNAlqRAWuiQVwkKXpEJY6JJUCGe5SFpTmpxZNsgMmlGYneMWuiQVwkKXpEIMVOgRcV5E3BIRt0XExXWFkiStXN+FHhHHAB8Cng+cDlwQEafXFUyStDKDbKGfBdyWmd/LzJ8D08CL64klSVqpyMz+7hjxMuC8zHxtdfuVwG9l5psWrbcN2FbdfDJwy6KHOgm4q68Qq8+szTBrM8xav2Hl/PXMPLnbSoNMW4wlxh7yv0Nm7gB2dHyQiNnMnBggx6oxazPM2gyz1m/Ucw6yy+UO4NS2248H7hwsjiSpX4MU+n8Dp0XEEyLi4cD5wNX1xJIkrVTfu1wy80hEvAn4CnAM8PHM3NPHQ3XcHTOCzNoMszbDrPUb6Zx9HxSVJI0W3ykqSYWw0CWpEI0VekTsi4hdEXFjRMxWYydGxLURcWv19YRqPCLib6tTCNwcEU9ve5wLq/VvjYgLG8i5ISI+GxHfiYi9EXHOiOZ8cvVczl9+GhFvHcWs1ff4o4jYExG7I+KKiHhkdQD969X3/XR1MJ2IeER1+7Zq+ea2x3l7NX5LRDyvoaxvqXLuiYi3VmMj8bxGxMcj4lBE7G4bqy1bRGytfk9vq+671HTkQbK+vHpefxkRE4vWX/JnGx1OKdLp9VNj1vdGqwdujogrI2LDKGRdkcxs5ALsA05aNPY3wMXV9YuB91TXXwB8idbc9rOBr1fjJwLfq76eUF0/oeaclwGvra4/HNgwijkXZT4G+BHw66OYFdgEfB9YV93+DHBR9fX8auzDwOur628APlxdPx/4dHX9dOAm4BHAE4DbgWNqznoGsBt4FK1JAv8KnDYqzyvwLODpwO4mfo+AbwDnVPf5EvD8mrM+hdYbCmeAibbxJX+21eV24Im0fh9vAk5vex095PVTY9ZzgWOr6+9pe16HmnVF/67GHnjpQr8F2Fhd3wjcUl3/CHDB4vWAC4CPtI0vWK+GjI+hVTwxyjmXyH0u8B+jmpVWof+QVoEcC1wDPI/WO+zmf2HOAb5SXf8KcE51/dhqvQDeDry97XEfXK/GrC8HPtZ2+8+BPx2l5xXYzMLiqSVbtew7beML1qsja9v4DAsLfcmfbfvron296vWw5Oun7qzVspcCl49K1l4vTe5DT+BfImJntN7+DzCemQcAqq+nVOPzBTDvjmqs03hdngj8GPhERHwrIj4WEcePYM7FzgeuqK6PXNbM3A+8D/gBcAA4DOwEfpKZR5b4vg9mqpYfBh67GllpbZ0/KyIeGxGPorWVeyoj+Ly2qSvbpur6amRebKVZH0vn108TXk3rL5a1kPVBTRb6MzPz6bTOxvjGiHjWMut2Oo1AT6cXGMCxtP7s+ofMfBpwH60/YTsZVs5fBWjti3sR8E/dVl1ibFWyVvt0X0zrz9PHAcfTeh10+r5Dy5qZe2n9eX0t8GVafzYfWeYuQ38NLGOl2YaZeWSzRsQ7aL0GLp8fWmGmoT2vjRV6Zt5ZfT0EXEnr7IwHI2IjQPX1ULV6p9MINH16gTuAOzLz69Xtz9Iq+FHL2e75wDcz82B1exSzPhf4fmb+ODN/AXweeAawISLm38zW/n0fzFQtXw/cvUpZycxLM/Ppmfms6vveymg+r/PqynZHdX01Mi+20qx30fn1U5vqgPHvAK/Ian/JqGZdSiOFHhHHR8Sj56/T2ue7m9apAeaPsF8IXFVdvxp4VXWU/mzgcPWn5FeAcyPihGqr79xqrBaZ+SPghxHx5GroOcC3Ry3nIhfwq90t85lGLesPgLMj4lHVrIn55/V64GUdss7/G14GfLX6ZboaOD9as2CeQOtg5TdqzkpEnFJ9/TXgd2k9v6P4vM6rJVu17N6IOLv6Ob2q7bGa1ulnu+QpRarXQ6fXTy0i4jzgbcCLMvP/RjlrR03smKe1b/qm6rIHeEc1/ljgOlpbQNcBJ1bjQevDMm4HdrHw4Mmrgduqyx80kPVMYBa4GfhnWrMARi5n9T0eBfwvsL5tbFSzvgv4Dq3/yD9Fa4bAE2n9ItxGa5fRI6p1H1ndvq1a/sS2x3lH9W+4hQFmYHTJ+u+0/sO5CXjOKD2vtP5zOQD8gtYW4WvqzAZMVD+j24G/Z9EEgRqyvrS6fj9wkIUHEZf82dI6jvHdatk72saXfP3UmPU2WvvEb6wuHx6FrCu5+NZ/SSqE7xSVpEJY6JJUCAtdkgphoUtSISx0SSqEhS5JhbDQJakQ/w+ABvFUWz0dmAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "point_estimates=[]\n",
    "for x in range(500):      #500次循环\n",
    "    sample=np.random.choice(a=data,size=100)    #每次随机抽样100个样本\n",
    "    point_estimates.append(sample.mean())\n",
    "    \n",
    "pd.DataFrame(point_estimates).hist(bins=40)    #均值大致呈钟型分布\n",
    "print('样本均值的均值为: ', np.array(point_estimates).mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 均值置信区间估计\n",
    "对于置信水平α与与之相对应的置信区间，置信区间中有α机会包含总体参数（如均值）。     \n",
    "这里使用课件中提供的定义函数并略加修改，可以指定相应的参数。    \n",
    "对于特定统计区间的人均GDP的均值，返回具有95%置信水平的置信区间，如下所示："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def makeConfidenceInterval(data, sample_size, _alpha):\n",
    "    sample=np.random.choice(a=data,size=sample_size)\n",
    "    sigma=sample.std()/(sample_size**0.5)\n",
    "    return stats.t.interval(alpha=_alpha,    #置信水平confidence level\n",
    "                 df=sample_size-1,  #自由度Degrees of freedom\n",
    "                loc=sample.mean(),\n",
    "                scale=sigma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(5848.845111513814, 11374.209888486186)\n"
     ]
    }
   ],
   "source": [
    "print(makeConfidenceInterval(data, 80, 0.95))\n",
    "#返回拥有95%置信水平的置信区间"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "修改置信度，查看不同置信水平下的置信区间大小："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "confidence 0.500000 has a interval of confidence range 2225.210000\n",
      "confidence 0.800000 has a interval of confidence range 3713.140000\n",
      "confidence 0.850000 has a interval of confidence range 3674.120000\n",
      "confidence 0.900000 has a interval of confidence range 5394.030000\n",
      "confidence 0.950000 has a interval of confidence range 5975.320000\n",
      "confidence 0.990000 has a interval of confidence range 6458.220000\n"
     ]
    }
   ],
   "source": [
    "for confidence in (.5, .8, .85, .9, .95, .99):\n",
    "    interval=makeConfidenceInterval(data, 80, confidence)\n",
    "    range_of_interval=round(interval[1]-interval[0],2)\n",
    "    print('confidence %f has a interval of confidence range %f' %(confidence,range_of_interval))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "之前计算出的均值的确包含在这个置信区间内的，不过也不是所有的时间内都会包含。    \n",
    "在重复运行中，出现过几次比较极端的区间，完全不包含实际的均值。     \n",
    "如果抽取的样本中有极端数据（比如美国），就有可能会导致区间不准确。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 置信区间验证\n",
    "由于置信区间不一定会包含要求的值，因此接下来，对于置信空间的准确性进行检验验证：     \n",
    "首先取出数据集的实际均值，然后多次计算置信区间，记录置信区间包含实际均值的频率。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9389\n"
     ]
    }
   ],
   "source": [
    "pop_mean=data.mean()\n",
    "times_in_interval=0;\n",
    "times = 10000\n",
    "for i in range(times):\n",
    "    interval=makeConfidenceInterval(data, 80, 0.95)\n",
    "    if (pop_mean>=interval[0]) & (pop_mean<=interval[1]):\n",
    "        times_in_interval += 1\n",
    "        \n",
    "print(times_in_interval/times)   #置信区间有95的可能性包含总体均值"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "改变置信水平，重新计算："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7977\n"
     ]
    }
   ],
   "source": [
    "pop_mean=data.mean()\n",
    "times_in_interval=0;\n",
    "times = 10000\n",
    "for i in range(times):\n",
    "    interval=makeConfidenceInterval(data, 80, 0.80)\n",
    "    if (pop_mean>=interval[0]) & (pop_mean<=interval[1]):\n",
    "        times_in_interval += 1\n",
    "        \n",
    "print(times_in_interval/times)   #置信区间有95的可能性包含总体均值"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "可以看出，置信区间的可信度与置信水平一致，这也正是“置信水平”一词的含义。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 随堂练习：Titanic-1\n",
    "对于不同舱位等级的分组，做个船费均值的区间估计："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "titanic = pd.read_csv(\"Titanic.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "自定义一个切分函数，这样的函数已经多次用到过了，此函数可以将指定的Pclass的旅客的所有信息放到一张新的DataFrame里面："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def slice_by_Pclass(dataframe, Pclass):\n",
    "    dataframe = dataframe[dataframe.Pclass==Pclass] \n",
    "    dataframe.columns.name = Pclass\n",
    "    return dataframe"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "将一等舱的乘客信息切分出来，看看切出来的数据如何："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'pandas.core.frame.DataFrame'>\n",
      "Int64Index: 216 entries, 1 to 889\n",
      "Data columns (total 12 columns):\n",
      "PassengerId    216 non-null int64\n",
      "Survived       216 non-null int64\n",
      "Pclass         216 non-null int64\n",
      "Name           216 non-null object\n",
      "Sex            216 non-null object\n",
      "Age            186 non-null float64\n",
      "SibSp          216 non-null int64\n",
      "Parch          216 non-null int64\n",
      "Ticket         216 non-null object\n",
      "Fare           216 non-null float64\n",
      "Cabin          176 non-null object\n",
      "Embarked       214 non-null object\n",
      "dtypes: float64(2), int64(5), object(5)\n",
      "memory usage: 17.7+ KB\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>1</th>\n",
       "      <th>PassengerId</th>\n",
       "      <th>Survived</th>\n",
       "      <th>Pclass</th>\n",
       "      <th>Name</th>\n",
       "      <th>Sex</th>\n",
       "      <th>Age</th>\n",
       "      <th>SibSp</th>\n",
       "      <th>Parch</th>\n",
       "      <th>Ticket</th>\n",
       "      <th>Fare</th>\n",
       "      <th>Cabin</th>\n",
       "      <th>Embarked</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>Cumings, Mrs. John Bradley (Florence Briggs Th...</td>\n",
       "      <td>female</td>\n",
       "      <td>38.0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>PC 17599</td>\n",
       "      <td>71.2833</td>\n",
       "      <td>C85</td>\n",
       "      <td>C</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>3</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>Futrelle, Mrs. Jacques Heath (Lily May Peel)</td>\n",
       "      <td>female</td>\n",
       "      <td>35.0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>113803</td>\n",
       "      <td>53.1000</td>\n",
       "      <td>C123</td>\n",
       "      <td>S</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>6</td>\n",
       "      <td>7</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>McCarthy, Mr. Timothy J</td>\n",
       "      <td>male</td>\n",
       "      <td>54.0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>17463</td>\n",
       "      <td>51.8625</td>\n",
       "      <td>E46</td>\n",
       "      <td>S</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>11</td>\n",
       "      <td>12</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>Bonnell, Miss. Elizabeth</td>\n",
       "      <td>female</td>\n",
       "      <td>58.0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>113783</td>\n",
       "      <td>26.5500</td>\n",
       "      <td>C103</td>\n",
       "      <td>S</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>23</td>\n",
       "      <td>24</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>Sloper, Mr. William Thompson</td>\n",
       "      <td>male</td>\n",
       "      <td>28.0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>113788</td>\n",
       "      <td>35.5000</td>\n",
       "      <td>A6</td>\n",
       "      <td>S</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "1   PassengerId  Survived  Pclass  \\\n",
       "1             2         1       1   \n",
       "3             4         1       1   \n",
       "6             7         0       1   \n",
       "11           12         1       1   \n",
       "23           24         1       1   \n",
       "\n",
       "1                                                Name     Sex   Age  SibSp  \\\n",
       "1   Cumings, Mrs. John Bradley (Florence Briggs Th...  female  38.0      1   \n",
       "3        Futrelle, Mrs. Jacques Heath (Lily May Peel)  female  35.0      1   \n",
       "6                             McCarthy, Mr. Timothy J    male  54.0      0   \n",
       "11                           Bonnell, Miss. Elizabeth  female  58.0      0   \n",
       "23                       Sloper, Mr. William Thompson    male  28.0      0   \n",
       "\n",
       "1   Parch    Ticket     Fare Cabin Embarked  \n",
       "1       0  PC 17599  71.2833   C85        C  \n",
       "3       0    113803  53.1000  C123        S  \n",
       "6       0     17463  51.8625   E46        S  \n",
       "11      0    113783  26.5500  C103        S  \n",
       "23      0    113788  35.5000    A6        S  "
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "new_data = slice_by_Pclass(titanic,1)\n",
    "new_data.info()\n",
    "new_data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接下来抽样，每次抽取100个样本，循环500次，记录每次的均值，然后输出："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "样本均值的均值为:  84.159923482\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAEOhJREFUeJzt3W+MXFd5x/HvAwFhspA/BEbGhJoKF0FtxcWrNAIJ7fKnCknVhApUUkqdErqoLS0Uv7F4A4giBZUQqS2q6ioRpqJZKCRNBJQqsrKkSIV2DQlrY0EomIAT2Y3iLNkQAU6fvphrZdns7NydP3tnz3w/0mhm7r1z5/HRnZ/Pnjn3TmQmkqTN72lNFyBJGgwDXZIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBXCQNdYi4gLI+K2iHgsIn4YEb/fdE1Sr85pugCpYZ8Afg60gN3AFyPi3sw82mxZ0vqFZ4pqXEXEucBpYGdmfrda9k/Aiczc32hxUg8cctE4+zXgibNhXrkX+PWG6pH6YqBrnE0AiyuWLQLPaaAWqW8GusbZEvDcFcueCzzaQC1S3wx0jbPvAudExI5lyy4B/EJUm5JfimqsRcQskMA7ac9y+RLwKme5aDOyh65x96fAFuAUcAvwJ4a5Nit76JJUCHvoklQIA12SCmGgS1IhDHRJKsSGXpzroosuyu3bt2/kW/bkscce49xzz226jJFnO9VjO9VjO3V2+PDhhzLz+d2229BA3759O/Pz8xv5lj2Zm5tjamqq6TJGnu1Uj+1Uj+3UWUT8sM52DrlIUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklQIA12SCmGgS1IhNvRMUY2H7fu/2HHd8euvHNl9S5udPXRJKoSBLkmF6BroEfGsiPiviLg3Io5GxIeq5S+JiK9HxH0R8ZmIeObwy5UkdVKnh/4z4LWZeQntX0W/PCIuAz4K3JiZO4DTwHXDK1OS1E3XQM+2perpM6pbAq8FPlctPwhcPZQKJUm1RGZ23yji6cBh4KXAJ4C/Br6WmS+t1l8M/Ftm7lzltTPADECr1dozOzs7uOqHZGlpiYmJiabLGHmd2mnhxGLH1+zadl5f7znMfQ+Lx1M9tlNn09PThzNzstt2taYtZuYTwO6IOB+4DXj5apt1eO0B4ADA5ORkboYL2Huh/Xo6tdO1a00tfNtTt1+PYe57WDye6rGd+reuWS6Z+QgwB1wGnB8RZ/9DeBHwwGBLkyStR51ZLs+veuZExBbg9cAx4C7gzdVme4Hbh1WkJKm7OkMuW4GD1Tj604DPZuYXIuLbwGxE/BXwTeCmIdYpSeqia6Bn5reA31hl+feBS4dRlCRp/TxTVJIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklQIA12SCmGgS1IhDHRJKoSBLkmFMNAlqRAGuiQVwkCXpEIY6JJUCANdkgphoEtSIQx0SSqEgS5JhTDQJakQBrokFcJAl6RCGOiSVAgDXZIK0TXQI+LiiLgrIo5FxNGIeE+1/IMRcSIi7qluVwy/XElSJ+fU2OYMsC8zvxERzwEOR8Sd1bobM/NjwytPklRX10DPzAeBB6vHj0bEMWDbsAuTJK1PZGb9jSO2A3cDO4H3AdcCPwHmaffiT6/ymhlgBqDVau2ZnZ3tt+ahW1paYmJioukyulo4sdhx3a5t5w39/Tu10zDrWmvf3WxEm6xmsxxPTbOdOpuenj6cmZPdtqsd6BExAXwF+Ehm3hoRLeAhIIEPA1sz8x1r7WNycjLn5+drvV+T5ubmmJqaarqMrrbv/2LHdcevv3Lo79+pnYZZ11r77mYj2mQ1m+V4aprt1FlE1Ar0WrNcIuIZwOeBT2fmrQCZeTIzn8jM/wP+Ebi0n4IlSf2pM8slgJuAY5n58WXLty7b7E3AkcGXJ0mqq84sl1cDbwcWIuKeatn7gWsiYjftIZfjwLuGUqEkqZY6s1y+CsQqq740+HIkSb2q00OXBqbbl5pNfXE5TAsnFrm24S+wNR489V+SCmGgS1IhDHRJKoSBLkmFMNAlqRAGuiQVwmmLY6rp68A0YRz/zRov9tAlqRAGuiQVwkCXpEIY6JJUCANdkgphoEtSIZy2OML6+bm1YTpb175dZ9a8iqCkjWUPXZIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBXCaYsSo/vj1aNal0aTPXRJKoSBLkmF6BroEXFxRNwVEcci4mhEvKdafmFE3BkR91X3Fwy/XElSJ3V66GeAfZn5cuAy4M8i4hXAfuBQZu4ADlXPJUkN6RromflgZn6jevwocAzYBlwFHKw2OwhcPawiJUndRWbW3zhiO3A3sBO4PzPPX7budGY+ZdglImaAGYBWq7Vndna2z5KHb2lpiYmJiabLYOHEYs+v3bXtvJ73Xfe1rS1w8vH117aWfuoepm51reXUw4trttNa++727+2nrlEzKp+7UTQ9PX04Mye7bVc70CNiAvgK8JHMvDUiHqkT6MtNTk7m/Px8rfdr0tzcHFNTU02X0dfVFrtNZ+vnB5OXX23xhoXBznztp+5h6md64N9++vY122mtfY/TtMVR+dyNooioFei1ZrlExDOAzwOfzsxbq8UnI2JrtX4rcKrXYiVJ/aszyyWAm4BjmfnxZavuAPZWj/cCtw++PElSXXX+Xn418HZgISLuqZa9H7ge+GxEXAfcD7xlOCVKkuroGuiZ+VUgOqx+3WDLkST1yjNFJakQXpxLqmGcZpto87KHLkmFMNAlqRAGuiQVwkCXpEIY6JJUCANdkgrhtEU9RVMXwJLUH3voklQIA12SCmGgS1IhDHRJKoSBLkmFcJaLRso4zrAZx3+zhsMeuiQVwkCXpEIY6JJUCANdkgphoEtSIQx0SSqE0xalAVhr6uG+XRtYiMaaPXRJKoSBLkmF6BroEXFzRJyKiCPLln0wIk5ExD3V7YrhlilJ6qZOD/2TwOWrLL8xM3dXty8NtixJ0np1DfTMvBt4eANqkST1oZ8x9HdHxLeqIZkLBlaRJKknkZndN4rYDnwhM3dWz1vAQ0ACHwa2ZuY7Orx2BpgBaLVae2ZnZwdS+DAtLS0xMTHRdBksnFhsuoQ1tbbAycebrmL0DbOddm07bzg7bsCofO5G0fT09OHMnOy2XU+BXnfdSpOTkzk/P9/1/Zo2NzfH1NRU02WM/GVV9+06ww0LnsrQzTDb6fj1Vw5lv00Ylc/dKIqIWoHe05BLRGxd9vRNwJFO20qSNkbXbkNE3AJMARdFxI+BDwBTEbGb9pDLceBdQ6xRklRD10DPzGtWWXzTEGqRJPXBM0UlqRAGuiQVwkCXpEIY6JJUCANdkgphoEtSIQx0SSqEgS5JhfBCHEM26tdjkVQOe+iSVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklQIA12SCmGgS1IhDHRJKoSBLkmFMNAlqRAGuiQVwkCXpEJ0DfSIuDkiTkXEkWXLLoyIOyPivur+guGWKUnqpk4P/ZPA5SuW7QcOZeYO4FD1XJLUoK6Bnpl3Aw+vWHwVcLB6fBC4esB1SZLWKTKz+0YR24EvZObO6vkjmXn+svWnM3PVYZeImAFmAFqt1p7Z2dkBlD1cS0tLTExMDGRfCycWB7KfUdTaAicfb7qK0ddkO+3adl4zb9yDQX7uSjM9PX04Mye7bTf0H4nOzAPAAYDJycmcmpoa9lv2bW5ujkHVeW3BPxK9b9cZbljwd8a7abKdjr9tqpH37cUgP3fjqtdZLicjYitAdX9qcCVJknrRa6DfAeytHu8Fbh9MOZKkXtWZtngL8J/AyyLixxFxHXA98IaIuA94Q/VcktSgrgN7mXlNh1WvG3AtkqQ+eKaoJBXCKQrA9hUzUfbtOvNLs1OOX3/lRpckNW7l52I9/Mw0wx66JBXCQJekQhjoklQIA12SCmGgS1IhDHRJKoSBLkmFMNAlqRAGuiQVwkCXpEIY6JJUCANdkgrhxbmkgq11gS0voFUee+iSVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpEE5blDRw3X6P1CmTw2EPXZIKYaBLUiH6GnKJiOPAo8ATwJnMnBxEUZKk9RvEGPp0Zj40gP1IkvrgkIskFSIys/cXR/wAOA0k8A+ZeWCVbWaAGYBWq7Vndna25/cbloUTi7/0vLUFTj7+5PNd286r/dpxsrKdtLpRbae1jmvY+GP7bDt1q2scTU9PH64zpN1voL8wMx+IiBcAdwJ/npl3d9p+cnIy5+fne36/YVk5xWrfrjPcsPDkaNRaU6y6Tc8q2cp20upGtZ26TR3c6GP7bDs5pfGpIqJWoPc15JKZD1T3p4DbgEv72Z8kqXc9B3pEnBsRzzn7GPgt4MigCpMkrU8/fwe2gNsi4ux+/jkzvzyQqiRJ69ZzoGfm94FLBliLJKkPTluUpEKM3lfvHXixH0lamz10SSqEgS5JhTDQJakQBrokFcJAl6RCGOiSVIhNM21RkmDtKczjPn3ZHrokFcJAl6RCGOiSVAgDXZIKYaBLUiHGYpZLvz+lNc4/M6dylXhc93sRv80+g8YeuiQVwkCXpEIY6JJUCANdkgphoEtSIQx0SSpEMdMWS5yCJY2jzfpZHoXfPbaHLkmFMNAlqRAGuiQVoq9Aj4jLI+I7EfG9iNg/qKIkSevXc6BHxNOBTwBvBF4BXBMRrxhUYZKk9emnh34p8L3M/H5m/hyYBa4aTFmSpPWKzOzthRFvBi7PzHdWz98O/GZmvnvFdjPATPX0ZcB3ei93w1wEPNR0EZuA7VSP7VSP7dTZr2Tm87tt1M889Fhl2VP+d8jMA8CBPt5nw0XEfGZONl3HqLOd6rGd6rGd+tfPkMuPgYuXPX8R8EB/5UiSetVPoP83sCMiXhIRzwTeCtwxmLIkSevV85BLZp6JiHcD/w48Hbg5M48OrLJmbaohogbZTvXYTvXYTn3q+UtRSdJo8UxRSSqEgS5JhRjrQI+Il0XEPctuP4mI90bEhRFxZ0TcV91f0HStTYuIv4yIoxFxJCJuiYhnVV+If71qp89UX46PtYh4T9VGRyPivdUyjycgIm6OiFMRcWTZslXbJtr+prqsyLci4pXNVb55jHWgZ+Z3MnN3Zu4G9gA/BW4D9gOHMnMHcKh6PrYiYhvwF8BkZu6k/SX4W4GPAjdW7XQauK65KpsXETuBP6Z9FvUlwG9HxA48ns76JHD5imWd2uaNwI7qNgP8/QbVuKmNdaCv8DrgfzLzh7QvYXCwWn4QuLqxqkbHOcCWiDgHeDbwIPBa4HPVetsJXg58LTN/mplngK8Ab8LjCYDMvBt4eMXiTm1zFfCpbPsacH5EbN2YSjcvA/1JbwVuqR63MvNBgOr+BY1VNQIy8wTwMeB+2kG+CBwGHqmCC9onmm1rpsKRcQR4TUQ8LyKeDVxB++Q7j6fOOrXNNuBHy7bz+KrBQAeqsd/fAf6l6VpGUTWueRXwEuCFwLm0/yReaaznwGbmMdrDUHcCXwbuBc6s+SJ1UuvSIvplBnrbG4FvZObJ6vnJs3/eVfenGqtsNLwe+EFm/m9m/gK4FXgV7T+Dz56c5qUfgMy8KTNfmZmvoT28cB8eT2vp1DZeWqQHBnrbNTw53ALtSxjsrR7vBW7f8IpGy/3AZRHx7IgI2t83fBu4C3hztY3tBETEC6r7FwO/S/u48njqrFPb3AH8YTXb5TJg8ezQjDob+zNFq7HOHwG/mpmL1bLnAZ8FXkw7zN6SmSu/zBkrEfEh4PdoDyF8E3gn7THNWeDCatkfZObPGityBETEfwDPA34BvC8zD3k8tUXELcAU7cvkngQ+APwrq7RN1XH4O9qzYn4K/FFmzjdR92Yy9oEuSaVwyEWSCmGgS1IhDHRJKoSBLkmFMNAlqRAGuiQVwkCXpEL8Pz8G0oMeibYBAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "point_estimates=[]\n",
    "\n",
    "for x in range(500):      #500次循环\n",
    "    sample=np.random.choice(a=new_data['Fare'],size=100)    #每次随机抽样100个样本\n",
    "    point_estimates.append(sample.mean())\n",
    "    \n",
    "pd.DataFrame(point_estimates).hist(bins=40)    #均值大致呈钟型分布\n",
    "print('样本均值的均值为: ', np.array(point_estimates).mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "最后，对所有三个舱位，输出均值的置信区间："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Class 1:  (71.54331363278226, 98.20568436721774)\n",
      "Class 2:  (16.158391739053343, 20.776190260946656)\n",
      "Class 3:  (12.410281464159652, 17.662126535840354)\n"
     ]
    }
   ],
   "source": [
    "interval_1 = makeConfidenceInterval(slice_by_Pclass(titanic,1)['Fare'], 100, 0.95)\n",
    "interval_2 = makeConfidenceInterval(slice_by_Pclass(titanic,2)['Fare'], 100, 0.95)\n",
    "interval_3 = makeConfidenceInterval(slice_by_Pclass(titanic,3)['Fare'], 100, 0.95)\n",
    "print('Class 1: ', interval_1)\n",
    "print('Class 2: ', interval_2)\n",
    "print('Class 3: ', interval_3)\n",
    "#返回拥有95%置信水平的置信区间"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 方差&标准差置信区间估计\n",
    "首先，随机抽取几组样本，将抽取的样本的方差与标准差作图，可以看出，值满足钟形分布："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[<matplotlib.axes._subplots.AxesSubplot object at 0x6baa5350>]],\n",
       "      dtype=object)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEVCAYAAADwyx6sAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAFSFJREFUeJzt3X2MZXV9x/H3V1iFMshDF6fbLbpaqRXZAu5ErTbNDGpLMRZMNZEaw1bM1rZam2xtiGmUak1oWkrTVFvXSkCjTK1KtQgqEabUWrCzCOzStaK4tSx2VwRWhhKbxW//uGfsdZjZe+7zPb99v5LJ3HvOuXM/c/bsZ86ch99EZiJJar4njTuAJGkwLHRJKoSFLkmFsNAlqRAWuiQVwkKXpEJY6JJUCAtdR7SIODkiro2IRyPiPyPi18edSerV0eMOII3Ze4H/BaaBs4DPRMSdmXn3eGNJ3QvvFNWRKiKOAx4CzsjMr1XTPgzsy8xLxhpO6oGHXHQk+xng8eUyr9wJPG9MeaS+WOg6kk0BB1dMOwgcP4YsUt8sdB3JloCnrpj2VOCRMWSR+mah60j2NeDoiDitbdqZgCdE1UieFNURLSLmgQTeSOsql+uBF3uVi5rIPXQd6X4bOBY4AFwD/JZlrqZyD12SCuEeuiQVwkKXpEJY6JJUCAtdkgox0sG51q9fn5s2beLRRx/luOOOG+VbD0QTczcxM5h71JqYu4mZobfcO3fufCAzT+m4YGYe9gM4BvgyrTEu7gb+qJp+FfBN4I7q46xOX2vLli2ZmXnzzTdnEzUxdxMzZ5p71JqYu4mZM3vLDSxmh37NzFp76N8HzsnMpYhYB3wxIm6o5r0tMz/e1Y8aSdJQdCz06qfDUvV0XfXhxeuSNGFqnRSNiKMi4g5ad9PdmJm3VbPeExF3RcQVEfGUoaWUJHXU1Z2iEXEicC3wFuC7wH8DTwZ2AN/IzHet8pptwDaA6enpLfPz8ywtLTE1NTWA+KPVxNxNzAzmHrUm5m5iZugt99zc3M7MnOm4YJ0D7e0fwDuB318xbRa4rtNrPSk6ek3MnGnuUWti7iZmzhzuSdGOh1wi4pRqz5yIOBZ4GfDViNhQTQvgAmB3Vz9yJEkDVecqlw3A1RFxFK1j7h/LzOsi4qaIOAUIWpctvmmIOSVJHdS5yuUu4OxVpp8zlESSpJ54678kFWKkt/6rDJsu+cxh5++97BUjSiKpnXvoklQIC12SCmGhS1IhLHRJKoSFLkmF8CoXFcOrb3Skcw9dkgphoUtSISx0SSqEhS5JhbDQJakQFrokFcJCl6RCWOiSVAgLXZIKYaFLUiEsdEkqhGO5SDgOjMrgHrokFcJCl6RCdCz0iDgmIr4cEXdGxN0R8UfV9GdGxG0RcU9E/F1EPHn4cSVJa6mzh/594JzMPBM4Czg3Il4E/AlwRWaeBjwEXDy8mJKkTjoWerYsVU/XVR8JnAN8vJp+NXDBUBJKkmqJzOy8UMRRwE7g2cB7gT8Fbs3MZ1fzTwVuyMwzVnntNmAbwPT09Jb5+XmWlpaYmpoa3HcxIk3MPYzMu/YdPOz8zRtP6Ps9esndT65BfU9N3EagmbmbmBl6yz03N7czM2c6LVfrssXMfBw4KyJOBK4FnrvaYmu8dgewA2BmZiZnZ2dZWFhgdna2zltPlCbmHkbmrZ0u8Xtd/+/XS+5+cg3qe2riNgLNzN3EzDDc3F1d5ZKZDwMLwIuAEyNi+QfCTwH3DzaaJKkbda5yOaXaMycijgVeBuwBbgZeXS12EfCpYYWUJHVW55DLBuDq6jj6k4CPZeZ1EfHvwHxE/DHwFeCDQ8wpSeqgY6Fn5l3A2atMvxd4wTBCSZK651guGqlOY6Ys27750KonKvsZU6Xue0tN5a3/klQIC12SCmGhS1IhLHRJKoSFLkmF8CoXDVyJV5P4F43UBO6hS1IhLHRJKoSFLkmFsNAlqRAWuiQVwkKXpEJY6JJUCAtdkgphoUtSISx0SSqEhS5JhXAsFz3BJI/FMsnZpHFzD12SCmGhS1IhOhZ6RJwaETdHxJ6IuDsi3lpNvzQi9kXEHdXHecOPK0laS51j6IeA7Zl5e0QcD+yMiBureVdk5p8NL54kqa6OhZ6Z3wa+XT1+JCL2ABuHHUyS1J2ujqFHxCbgbOC2atKbI+KuiLgyIk4acDZJUhciM+stGDEF/BPwnsz8ZERMAw8ACbwb2JCZb1jldduAbQDT09Nb5ufnWVpaYmpqalDfw8g0MXcvmXftOzikNPVNHwv7Hxt3ivo2bzwBaOY2As3M3cTM0Fvuubm5nZk502m5WoUeEeuA64DPZeafrzJ/E3BdZp5xuK8zMzOTi4uLLCwsMDs72/F9J00Tc/eSeRKu9d6++RCX72rObRLLf1O0idsINDN3EzNDb7kjolah17nKJYAPAnvayzwiNrQt9ipgd1cJJUkDVWcX6CXA64FdEXFHNe3twIURcRatQy57gd8cSkJJUi11rnL5IhCrzLp+8HEkSb3yTlFJKoSFLkmFsNAlqRAWuiQVwkKXpEJY6JJUCAtdkgphoUtSISx0SSqEhS5JhbDQJakQFrokFcJCl6RCWOiSVAgLXZIKYaFLUiEsdEkqhIUuSYWw0CWpEHX+SLSkDjZd8hkAtm8+xNbq8bK9l71iHJF0BHIPXZIKYaFLUiE6FnpEnBoRN0fEnoi4OyLeWk0/OSJujIh7qs8nDT+uJGktdfbQDwHbM/O5wIuA34mI04FLgC9k5mnAF6rnkqQx6VjomfntzLy9evwIsAfYCJwPXF0tdjVwwbBCSpI6i8ysv3DEJuAW4AzgW5l5Ytu8hzLzCYddImIbsA1genp6y/z8PEtLS0xNTfUZffSamLuXzLv2HRxSmvqmj4X9j407RfdWy7154wlDe79O/1Z13/tI2bYnQS+55+bmdmbmTKflahd6REwB/wS8JzM/GREP1yn0djMzM7m4uMjCwgKzs7O13neSNDF3L5k3rbjsbhy2bz7E5buad1XtarmHedlip3+ruu99pGzbk6CX3BFRq9BrXeUSEeuATwAfycxPVpP3R8SGav4G4EBXCSVJA1XnKpcAPgjsycw/b5v1aeCi6vFFwKcGH0+SVFed32lfArwe2BURd1TT3g5cBnwsIi4GvgW8ZjgRJUl1dCz0zPwiEGvMfulg40iSetW8s05SwwzqxKXUibf+S1IhLHRJKoSFLkmFsNAlqRAWuiQVwkKXpEJY6JJUCAtdkgphoUtSISx0SSqEhS5JhXAsl0Itjx+yffMhtq4ylojjhxwZDjeOjNtAedxDl6RCWOiSVAgLXZIKYaFLUiEsdEkqhFe5SGPW6S8aSXW5hy5JhbDQJakQHQs9Iq6MiAMRsbtt2qURsS8i7qg+zhtuTElSJ3X20K8Czl1l+hWZeVb1cf1gY0mSutWx0DPzFuDBEWSRJPUhMrPzQhGbgOsy84zq+aXAVuB7wCKwPTMfWuO124BtANPT01vm5+dZWlpiampqAPFHq0m5d+07CMD0sbD/sSfO37zxhI6vHae1ck+6JuVu3waatG0va2Jm6C333Nzczsyc6bRcr4U+DTwAJPBuYENmvqHT15mZmcnFxUUWFhaYnZ3t+L6Tpkm52wfnunzXE69OPdzATJNwGd1auSddk3K3bwNN2raXNTEz9JY7ImoVek9XuWTm/sx8PDN/AHwAeEEvX0eSNDg9FXpEbGh7+ipg91rLSpJGo+PvhhFxDTALrI+I+4B3ArMRcRatQy57gd8cYkZJUg0dCz0zL1xl8geHkEWS1AfvFJWkQjTjdLwGbhKuZJE0WO6hS1IhLHRJKoSFLkmFsNAlqRAWuiQVwkKXpEJY6JJUCAtdkgphoUtSISx0SSqEhS5JhbDQJakQFrokFcJCl6RCWOiSVAgLXZIKYaFLUiEsdEkqhIUuSYWw0CWpEB0LPSKujIgDEbG7bdrJEXFjRNxTfT5puDElSZ3U2UO/Cjh3xbRLgC9k5mnAF6rnkqQx6ljomXkL8OCKyecDV1ePrwYuGHAuSVKXIjM7LxSxCbguM8+onj+cmSe2zX8oM1c97BIR24BtANPT01vm5+dZWlpiampqAPFHa9Jy79p3sOMy08fC/sdGEGbAzD18mzee8MPHk7Zt19HEzNBb7rm5uZ2ZOdNpuaN7TlVTZu4AdgDMzMzk7OwsCwsLzM7ODvutB27Scm+95DMdl9m++RCX7xr6P/PAmXv49r5u9oePJ23brqOJmWG4uXu9ymV/RGwAqD4fGFwkSVIvei30TwMXVY8vAj41mDiSpF7VuWzxGuBfgedExH0RcTFwGfDyiLgHeHn1XJI0Rh0P9mXmhWvMeumAs0iS+tCMszeSRm7TYU66773sFSNMorq89V+SCmGhS1IhLHRJKoSFLkmFsNAlqRBe5SKpa4e7Aga8CmZc3EOXpEJY6JJUCAtdkgphoUtSISx0SSqEV7mMkVcKaJzat7/tmw/V+oMpmmzuoUtSISx0SSqEhS5JhbDQJakQFrokFcJCl6RCWOiSVAgLXZIK0deNRRGxF3gEeBw4lJkzgwglSereIO4UncvMBwbwdSRJffCQiyQVIjKz9xdHfBN4CEjg/Zm5Y5VltgHbAKanp7fMz8+ztLTE1NRUz+87LoPOvWvfwcPO37zxhL5eDzB9LOx/rKtYE8HcozXo3J223UE4knpkbm5uZ51D2v0W+k9m5v0R8TTgRuAtmXnLWsvPzMzk4uIiCwsLzM7O9vy+4zLo3P0OztXp9dAadOnyXc0bg83cozXo3KMYWO5I6pGIqFXofR1yycz7q88HgGuBF/Tz9SRJveu50CPiuIg4fvkx8EvA7kEFkyR1p5/fsaaBayNi+et8NDM/O5BUkqSu9VzomXkvcOYAs0iS+uBli5JUCAtdkgphoUtSISx0SSqEhS5JhbDQJakQFrokFaJ5g06MwfKYKds3H2JrjfFT2vUzpkWdsVokaZl76JJUCAtdkgphoUtSISx0SSqEhS5JhfAqlyHzShWpO3X/zwzjqrPDvfco/gpTv9xDl6RCWOiSVAgLXZIKYaFLUiEsdEkqRGOucunnapF+zmxL6l6n/1PjumKkn//r/fbEKL5n99AlqRAWuiQVoq9Cj4hzI+I/IuLrEXHJoEJJkrrXc6FHxFHAe4FfAU4HLoyI0wcVTJLUnX720F8AfD0z783M/wXmgfMHE0uS1K3IzN5eGPFq4NzMfGP1/PXACzPzzSuW2wZsq54+B/gPYD3wQK+hx6iJuZuYGcw9ak3M3cTM0FvuZ2TmKZ0W6ueyxVhl2hN+OmTmDmDHj7wwYjEzZ/p477FoYu4mZgZzj1oTczcxMww3dz+HXO4DTm17/lPA/f3FkST1qp9C/zfgtIh4ZkQ8GXgt8OnBxJIkdavnQy6ZeSgi3gx8DjgKuDIz76758h2dF5lITczdxMxg7lFrYu4mZoYh5u75pKgkabJ4p6gkFcJCl6RCDLTQOw0FEBFbI+I7EXFH9fHGtnkXRcQ91cdFg8w1gNxXtGX+WkQ83Dbv8bZ5IzspHBFXRsSBiNi9xvyIiL+svqe7IuL5bfPGua475X5dlfeuiPhSRJzZNm9vROyq1vXi6FLXyj0bEQfbtoV3tM0byxAZNTK/rS3v7mpbPrmaN851fWpE3BwReyLi7oh46yrLTNz2XTP3cLfvzBzIB60To98AngU8GbgTOH3FMluBv1rltScD91afT6oenzSobP3mXrH8W2idAF5+vjSKnKvk+EXg+cDuNeafB9xA636BFwG3jXtd18z94uU8tIaVuK1t3l5g/YSu71ngun63r1FmXrHsK4GbJmRdbwCeXz0+HvjaKl0ycdt3zdxD3b4HuYfez1AAvwzcmJkPZuZDwI3AuQPMdjjd5r4QuGYkyQ4jM28BHjzMIucDH8qWW4ETI2ID413XHXNn5peqXAC30rq/YexqrO+1jG2IjC4zT8R2DZCZ387M26vHjwB7gI0rFpu47btO7mFv34Ms9I3Af7U9v48n/iMA/Fr168bHI2L5xqS6rx2G2u8dEc8Angnc1Db5mIhYjIhbI+KC4cXs2lrf1zjXdbcuprUXtiyBz0fEzmgNKTFpfj4i7oyIGyLiedW0iV/fEfFjtErvE22TJ2JdR8Qm4GzgthWzJnr7PkzudgPfvgf5F4vqDAXwj8A1mfn9iHgTcDVwTs3XDks37/1a4OOZ+XjbtKdn5v0R8SzgpojYlZnfGHjK7q31fY1zXdcWEXO0NvhfaJv8kmpdPw24MSK+Wu2FToLbaY23sRQR5wH/AJxGM9b3K4F/ycz2vfmxr+uImKL1Q+b3MvN7K2ev8pKJ2L475F5eZijb9yD30DsOBZCZ383M71dPPwBsqfvaIermvV/Lil9LM/P+6vO9wAKtn8qTYK3va+KHbIiInwP+Fjg/M7+7PL1tXR8ArqV1OGMiZOb3MnOpenw9sC4i1tOA9c3ht+uxrOuIWEerFD+SmZ9cZZGJ3L5r5B7u9j3AEwJH0zoB8Uz+/+TP81aeNGh7/Crg1vz/ExnfpHUS46Tq8cmDytZv7mq559A6aRFt004CnlI9Xg/cw4hOeFXvuYm1T9K9gh89afTlca/rmrmfDnwdePGK6ccBx7c9/hKt0T4nJfdPLG8b1X/Eb1Xrvtb2NY7M1fwTaB1nP25S1nW13j4E/MVhlpm47btm7qFu3wM75JJrDAUQEe8CFjPz08DvRsSvAoeqjWhr9doHI+LdtMaHAXhX/uivf0NTMze0ThrNZ7XGK88F3h8RP6D1285lmfnvo8gdEdfQurJifUTcB7wTWFd9T38DXE/rSoCvA/8D/EY1b2zrumbudwA/DrwvIgAOZWtkumng2mra0cBHM/OzE5T71cBvRcQh4DHgtdW20s8QGcPODK0dq89n5qNtLx3rugZeArwe2BURd1TT3k6rDCd5+66Te6jbt7f+S1IhvFNUkgphoUtSISx0SSqEhS5JhbDQJWlIOg2QtmLZp1eDe32lupv+vG7fz0KXpOG5ivpjyfwh8LHMPJvWzV7v6/bNLHRJGpJcZYC0iPjpiPhsNWbLP0fEzy4vDjy1enwCPdzhOsixXCRJne0A3pSZ90TEC2ntiZ8DXEprcK630Lpb9GXdfmELXZJGpBq468XA31d3hQI8pfp8IXBVZl4eET8PfDgizsjMH9T9+ha6JI3Ok4CHM/OsVeZdTHW8PTP/NSKOoTVG1IFuvrgkaQSyNZzuNyPiNfDDP6W3/GfovgW8tJr+XOAY4DvdfH3HcpGkIWkfIA3YT2uAtJuAv6b1J+vW0Rr0710RcTqtYcWnaJ0g/YPM/HxX72ehS1IZPOQiSYWw0CWpEBa6JBXCQpekQljoklQIC12SCmGhS1Ih/g+PkgsMtCbNnQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "point_estimates=[]\n",
    "for x in range(500):      #500次循环\n",
    "    sample=np.random.choice(a=data,size=100)    #每次随机抽样100个样本\n",
    "    point_estimates.append(sample.var())\n",
    "    \n",
    "pd.DataFrame(point_estimates).hist(bins=40)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[<matplotlib.axes._subplots.AxesSubplot object at 0x6ba60050>]],\n",
       "      dtype=object)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAFW9JREFUeJzt3X+M5Hd93/HnG2Pg4gX/wGZ7XJweBDfCeIXhVsRJJLRrIDGmio2aSFCE7EK0aQiURJe2l1hKTAmSCThISVCRI7s4CLw4BAtqIHB1vbhIAbpHz967HOBfF+LDOde1fWGJ5WTNu3/M93Jz65md2Znv7Mz3c8+HNNqZ7/c733ndfGdf993vr4nMRJLUfM8adwBJUj0sdEkqhIUuSYWw0CWpEBa6JBXCQpekQljoklQIC12ntIg4JyJui4gfRsTfRMS/HXcmaVDPHncAacw+CvwjMA1cDHwhIu7OzIPjjSVtXnimqE5VEXEG8DhwUWZ+txr2CeBIZu4ZazhpAG5y0ansXwFPHy/zyt3AK8aURxqKha5T2RRwbN2wY8Dzx5BFGpqFrlPZKvCCdcNeAPxgDFmkoVnoOpV9F3h2RFzQNuyVgDtE1UjuFNUpLSIWgQR+hdZRLl8EftajXNRErqHrVPcuYBvwCHAL8GuWuZrKNXRJKoRr6JJUCAtdkgphoUtSISx0SSrEll6c69xzz82dO3eObP4//OEPOeOMM0Y2/zo1Jas569eUrE3JCc3JOmjOffv2PZqZ5/WcMDO37LZr164cpTvvvHOk869TU7Kas35NydqUnJnNyTpoTmA5++hYN7lIUiEsdEkqhIUuSYWw0CWpEBa6JBXCQpekQljoklQIC12SCmGhS1IhtvTUf6mpdu75wobjD1/3pi1KInXnGrokFcJCl6RCWOiSVAgLXZIKYaFLUiEsdEkqhIUuSYWw0CWpEBa6JBXCQpekQljoklQIC12SCtGz0CPieRHxzYi4OyIORsT7quEfj4gHI2J/dbt49HElSd30c7XFp4BLM3M1Ik4HvhYRX6rG/cfM/Mzo4kmS+tWz0DMzgdXq4enVLUcZSpK0edHq6x4TRZwG7ANeBnw0M/9zRHwc+Blaa/B3AHsy86kOz10AFgCmp6d3LS4u1pd+ndXVVaampkY2/zo1Jas5W1aOHNtw/MyOM/uel+9p/ZqSddCc8/Pz+zJzttd0fRX6P08ccRZwG/Ae4P8Bfwc8B7gBuD8z/8tGz5+dnc3l5eW+X2+zlpaWmJubG9n869SUrOZsqfMLLnxP69eUrIPmjIi+Cn1TR7lk5hPAEnBZZj6cLU8B/w14zaZTSpJq089RLudVa+ZExDbg9cC3I2J7NSyAK4EDowwqSdpYP0e5bAdurrajPwu4NTNvj4j/GRHnAQHsB/79CHNKknro5yiXe4BXdRh+6UgSSZIG4pmiklQIC12SCmGhS1IhLHRJKoSFLkmFsNAlqRAWuiQVwkKXpEJY6JJUCAtdkgrRz7VcpEbodYlbqXSuoUtSISx0SSqEhS5JhbDQJakQFrokFcKjXKQabHSEzWa+QFoahmvoklSIfr4k+nkR8c2IuDsiDkbE+6rhL4mIb0TEvRHx6Yh4zujjSpK66WcN/Sng0sx8JXAxcFlEXAJ8EPhIZl4APA68c3QxJUm99Cz0bFmtHp5e3RK4FPhMNfxm4MqRJJQk9SUys/dEEacB+4CXAR8FPgR8PTNfVo0/H/hSZl7U4bkLwALA9PT0rsXFxfrSr7O6usrU1NTI5l+npmSdtJwrR451HD69DY4+ucVh+jSz48yTHk/ae9pNU3JCc7IOmnN+fn5fZs72mq6vo1wy82ng4og4C7gNeHmnybo89wbgBoDZ2dmcm5vr5yUHsrS0xCjnX6emZJ20nFd3OZpk98wa169M5kFbh982d9LjSXtPu2lKTmhO1lHn3NRRLpn5BLAEXAKcFRHHf4N+HPh+vdEkSZvRz1Eu51Vr5kTENuD1wCHgTuCXqsmuAj43qpCSpN76+Rt1O3BztR39WcCtmXl7RPw1sBgRvw/8H+DGEeaUJPXQs9Az8x7gVR2GPwC8ZhShJEmbN5l7kaSCrL8swO6ZtZN27nppANXFU/8lqRAWuiQVwkKXpEJY6JJUCAtdkgphoUtSISx0SSqEhS5JhbDQJakQFrokFcJCl6RCWOiSVAgLXZIKYaFLUiEsdEkqhIUuSYWw0CWpEP18SfT5EXFnRByKiIMR8d5q+LURcSQi9le3y0cfV5LUTT9fQbcG7M7Mb0XE84F9EbG3GveRzPzw6OJJkvrVz5dEPww8XN3/QUQcAnaMOpgkaXM2tQ09InYCrwK+UQ16d0TcExE3RcTZNWeTJG1CZGZ/E0ZMAV8FPpCZn42IaeBRIIH3A9sz8x0dnrcALABMT0/vWlxcrCv7M6yurjI1NTWy+depKVknLefKkWMdh09vg6NPbnGYAW0m68yOM0cbZgOTtuw30pSsg+acn5/fl5mzvabrq9Aj4nTgduDLmfmHHcbvBG7PzIs2ms/s7GwuLy/3fL1BLS0tMTc3N7L516kpWSct5849X+g4fPfMGtev9LNLaPw2k/XwdW8acZruJm3Zb6QpWQfNGRF9FXo/R7kEcCNwqL3MI2J722RvBg5sOqUkqTb9rCb8HPB2YCUi9lfDfgd4a0RcTGuTy2HgV0eSUJLUl36OcvkaEB1GfbH+OJKkQXmmqCQVwkKXpEJY6JJUCAtdkgphoUtSISx0SSqEhS5JhbDQJakQFrokFcJCl6RCWOiSVIhmXG9UOkV1u1zwceO8vK4mj2voklQIC12SCmGhS1IhLHRJKoSFLkmFsNAlqRAWuiQVomehR8T5EXFnRByKiIMR8d5q+DkRsTci7q1+nj36uJKkbvpZQ18Ddmfmy4FLgF+PiAuBPcAdmXkBcEf1WJI0Jj0LPTMfzsxvVfd/ABwCdgBXADdXk90MXDmqkJKk3iIz+584YidwF3AR8L3MPKtt3OOZ+YzNLhGxACwATE9P71pcXBwycnerq6tMTU2NbP51akrWScu5cuRYx+HT2+Dok1scZkBbmXVmx5kDP3fSlv1GmpJ10Jzz8/P7MnO213R9F3pETAFfBT6QmZ+NiCf6KfR2s7Ozuby83NfrDWJpaYm5ubmRzb9OTck6aTm7Xdtk98wa168049JEW5l1mGu9TNqy30hTsg6aMyL6KvS+jnKJiNOBvwA+mZmfrQYfjYjt1fjtwCObTilJqk0/R7kEcCNwKDP/sG3U54GrqvtXAZ+rP54kqV/9/N33c8DbgZWI2F8N+x3gOuDWiHgn8D3gl0cTUZLUj56FnplfA6LL6NfVG0eSNCjPFJWkQljoklQIC12SCmGhS1IhLHRJKkQzTq1TMfwW+6210fvte10e19AlqRAWuiQVwkKXpEJY6JJUCHeKatNGuWOz17wldecauiQVwkKXpEJY6JJUCAtdkgphoUtSITzKRVJHK0eOcXWXo468bMBkcg1dkgrRz5dE3xQRj0TEgbZh10bEkYjYX90uH21MSVIv/ayhfxy4rMPwj2TmxdXti/XGkiRtVs9Cz8y7gMe2IIskaQiRmb0nitgJ3J6ZF1WPrwWuBv4eWAZ2Z+bjXZ67ACwATE9P71pcXKwhdmerq6tMTU2NbP51akrWTjlXjhwbU5ruprfB0SfHnaI/k5J1ZseZG45/5LFjXXP2eu5Wa/LvUz/m5+f3ZeZsr+kGLfRp4FEggfcD2zPzHb3mMzs7m8vLyz1fb1BLS0vMzc2NbP51akrWTjkn8Xoru2fWuH6lGQdtTUrWXkeq/PEnP9c156Qd5dLk36d+RERfhT7QUS6ZeTQzn87MHwF/CrxmkPlIkuozUKFHxPa2h28GDnSbVpK0NXr+3RcRtwBzwLkR8RDwe8BcRFxMa5PLYeBXR5hRktSHnoWemW/tMPjGEWSRJA3BM0UlqRAWuiQVwkKXpEJY6JJUCAtdkgox/tPVNBYbne05aWcBSuqPa+iSVAgLXZIKYaFLUiEsdEkqhDtFpVNUr8sg754Z/LnuWB8P19AlqRAWuiQVwkKXpEJY6JJUCAtdkgrhUS56hvYjGHbPrHH1BH4ptKRncg1dkgrRs9Aj4qaIeCQiDrQNOyci9kbEvdXPs0cbU5LUSz9r6B8HLls3bA9wR2ZeANxRPZYkjVHPQs/Mu4DH1g2+Ari5un8zcGXNuSRJmxSZ2XuiiJ3A7Zl5UfX4icw8q23845nZcbNLRCwACwDT09O7FhcXa4jd2erqKlNTUyObf53GnXXlyLG+ppveBkefHHGYGjQlJzQn6yhzzuw4s9b5jfv3qV+D5pyfn9+XmbO9phv5US6ZeQNwA8Ds7GzOzc2N7LWWlpYY5fzrNO6s/R65sntmjetXJv9gqKbkhOZkHWXOw2+bq3V+4/596teocw56lMvRiNgOUP18pL5IkqRBDFronweuqu5fBXyunjiSpEH1c9jiLcBfAT8VEQ9FxDuB64A3RMS9wBuqx5KkMeq5gSwz39pl1OtqziJJGoJnikpSISx0SSqEhS5JhbDQJakQFrokFcJCl6RCWOiSVAgLXZIKYaFLUiEsdEkqhIUuSYWw0CWpEBa6JBXCQpekQljoklQIC12SCjH531R7CtvZ5xc5SxK4hi5JxRhqDT0iDgM/AJ4G1jJzto5QkqTNq2OTy3xmPlrDfCRJQ3CTiyQVYthCT+ArEbEvIhbqCCRJGkxk5uBPjnhxZn4/Il4E7AXek5l3rZtmAVgAmJ6e3rW4uDhM3g2trq4yNTU1svnXqZ+sK0eObVGa7qa3wdEnx52it6bkhOZkHWXOmR1ndh3X63Pf6blN+d0fNOf8/Py+fvZRDlXoJ80o4lpgNTM/3G2a2dnZXF5eruX1OllaWmJubm5k869TP1kn4bDF3TNrXL8y+Ue3NiUnNCfrKHMevu5NXcf1+tx3em5TfvcHzRkRfRX6wJtcIuKMiHj+8fvAzwMHBp2fJGk4w/z3Ow3cFhHH5/OpzPzLWlJJkjZt4ELPzAeAV9aYRZI0hMnfkCdJbTptY989s8bVe76w4bb5U4HHoUtSISx0SSqEhS5JhbDQJakQ7hTtw0YnOvTaCdPtubtn1pgbJpTUYKM6aW6Qk5JKeO3jXEOXpEJY6JJUCAtdkgphoUtSISx0SSqER7lIUmUSLlk9DNfQJakQFrokFcJCl6RCWOiSVIjG7BSdhNNqOxlmJ0rTd8BITVP675xr6JJUCAtdkgoxVKFHxGUR8Z2IuC8i9tQVSpK0eQMXekScBnwUeCNwIfDWiLiwrmCSpM0ZZg39NcB9mflAZv4jsAhcUU8sSdJmRWYO9sSIXwIuy8xfqR6/HfjpzHz3uukWgIXq4U8B3xk8bk/nAo+OcP51akpWc9avKVmbkhOak3XQnP8yM8/rNdEwhy1Gh2HP+N8hM28AbhjidfoWEcuZObsVrzWspmQ1Z/2akrUpOaE5WUedc5hNLg8B57c9/nHg+8PFkSQNaphC/9/ABRHxkoh4DvAW4PP1xJIkbdbAm1wycy0i3g18GTgNuCkzD9aWbDBbsmmnJk3Jas76NSVrU3JCc7KONOfAO0UlSZPFM0UlqRAWuiQVYuILPSJ+MyIORsSBiLglIp5X7Yj9RkTcGxGfrnbKEhHPrR7fV43f2Taf366GfycifmFEWd9b5TwYEb9RDTsnIvZWWfdGxNnV8IiIP6oy3RMRr26bz1XV9PdGxFU15LopIh6JiANtw2rLFRG7ImKles4fRUSnQ1qHyfrL1Xv6o4iYXTd9x+Xa7bIU3T47NeX8UER8u3rfbouIs8adc4Os769y7o+Ir0TEi6vhY1v+nXK2jfutiMiIOHcSc0bEtRFxpHo/90fE5W3jtm7ZZ+bE3oAdwIPAturxrcDV1c+3VMM+Bvxadf9dwMeq+28BPl3dvxC4G3gu8BLgfuC0mrNeBBwAfozWzub/AVwA/AGwp5pmD/DB6v7lwJdoHc9/CfCNavg5wAPVz7Or+2cPme21wKuBA23DassFfBP4meo5XwLeWHPWl9M6KW0JmG0b3nG5Vrf7gZcCz6mmubDtM/SMz05NOX8eeHZ1/4Nt7+nYcm6Q9QVt9/8DJ35vxrb8O+Wshp9P6+CLvwHOncScwLXAb3WYdkuX/cSvodMqx20R8WxaZfkwcCnwmWr8zcCV1f0rqsdU419X/S98BbCYmU9l5oPAfbQuXVCnlwNfz8x/yMw14KvAm9dlWp/1z7Ll68BZEbEd+AVgb2Y+lpmPA3uBy4YJlpl3AY+tG1xLrmrcCzLzr7L1CfyztnnVkjUzD2VmpzOMuy3XjpelqD4L3T47deT8SrXsAb5O69yMsebcIOvftz08gxMnBY5t+Xf5nAJ8BPhPnHzi4iTm7GRLl/1EF3pmHgE+DHyPVpEfA/YBT7T94jxEa02e6uffVs9dq6Z/YfvwDs+pywHgtRHxwoj4MVprEOcD05n5cJXpYeBF67Ouy7QVWakx147q/qjzdrLZrC+k+2enbu+gtRY4sTkj4gMR8bfA24DfHTDrSJd/RPwicCQz7143aqJyVt5dbf65KapNmAPkHGrZT3ShV2/KFbT+VHkxrTWJN3aY9Pj/3N0uR9DXZQqGkZmHaP2ZvRf4S1p/Qq1t8JSxZe1hs7nGmXcis0bENbSW/SePD9pkni3JmZnXZOb5tHIevwbTxGStVoyu4cR/NieN3mSeUb+n/xX4SeBiWiuf11fDtzTnRBc68Hrgwcz8v5n5T8BngZ+l9efV8ZOi2i858M+XI6jGn0nrT6MtuUxBZt6Yma/OzNdWr3svcLT6c4/q5yPrs67LtFWXVKgr10Oc2LQwyrydbDbro3T/7NSi2gn3r4G3VX/aT2TOdT4F/JsBs45y+f8krZW5uyPicDXvb0XEv5iwnGTm0cx8OjN/BPwpJzbpbu2yH2SnwFbdgJ8GDtLadh60tie9B/hzTt5p8K7q/q9z8k7RW6v7r+DkHRMPUPNO0ep1XlT9/Ang27R2ynyIk3c+/kF1/02cvFPnm3lip86D1XPPru6fU0O2nZy8E6e2XLQuA3EJJ3Y2XV5n1rbhS5y8U7TjcqW13+WBatjxHU6vqJ7T8bNT03t6GfDXwHnrphtrzi5ZL2i7/x7gM5Ow/Lst+2rcYU7sFJ2onMD2tvu/SWu7+ZYv+1oLbRQ34H20yvEA8InqjXkprT3W91X/+OdW0z6venxfNf6lbfO5htZe5e8wxFEYPbL+r+oX+m7gddWwFwJ30Fpbv6PtwxW0viDkfmCFk4vqHdW/4T7g39WQ6xZafwb+E601g3fWmQuYrZbP/cCfUJ2BXGPWN1f3nwKOAl/utVxp7cP4bjXumrbhHT87NeW8j9Z20f3V7WPjzrlB1r+oltk9wH8Hdox7+XfKuW78YU4U+kTlpNVNK9X7+XlOLvgtW/ae+i9JhZj0beiSpD5Z6JJUCAtdkgphoUtSISx0SSqEhS5JhbDQJakQ/x8LHZtQaabvTgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "point_estimates=[]\n",
    "for x in range(500):      #500次循环\n",
    "    sample=np.random.choice(a=data,size=100)    #每次随机抽样100个样本\n",
    "    point_estimates.append(sample.std())\n",
    "    \n",
    "pd.DataFrame(point_estimates).hist(bins=40)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "下面进行置信区间的计算。     \n",
    "对于置信水平α与与之相对应的置信区间，置信区间中有α机会包含总体参数（如均值）。    \n",
    "通过查阅资料可以得知，对于方差的置信区间，可以由以下公式求得：\n",
    "$$\\left ( \\frac{\\left ( n-1 \\right )S^{2}}{\\chi_{\\alpha /2}^{2}\\left ( n-1 \\right ) },\\frac{\\left ( n-1 \\right )S^{2}}{\\chi_{1-\\alpha /2}^{2}\\left ( n-1 \\right ) } \\right )$$\n",
    "而对于标准差，只需要将区间值开方即可。       \n",
    "对于特定统计区间的人均GDP的方差与标准差，返回具有95%置信水平的置信区间，如下所示："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "def makeConfidenceInterval_var(data, sample_size, _alpha):\n",
    "    _alpha = 1 - _alpha\n",
    "    sample=np.random.choice(a=data,size=sample_size)\n",
    "    sigma = sample.var() * (sample_size - 1)\n",
    "    chi1 = stats.chi2.ppf((1 - _alpha / 2), sample_size - 1)\n",
    "    chi2 = stats.chi2.ppf(_alpha / 2, sample_size - 1)\n",
    "    interval = [sigma / chi1, sigma / chi2]\n",
    "    return interval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[138866624.64604267, 260112416.51257148]\n"
     ]
    }
   ],
   "source": [
    "print(makeConfidenceInterval_var(data, 80, 0.95))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "def makeConfidenceInterval_std(data, sample_size, _alpha):\n",
    "    _alpha = 1 - _alpha\n",
    "    sample=np.random.choice(a=data,size=sample_size)\n",
    "    sigma = sample.std() * ((sample_size - 1) ** 0.5)\n",
    "    chi1 = stats.chi2.ppf((1 - _alpha / 2), sample_size - 1)\n",
    "    chi2 = stats.chi2.ppf(_alpha / 2, sample_size - 1)\n",
    "    interval = [sigma / (chi1 ** 0.5), sigma / (chi2 ** 0.5)]\n",
    "    return interval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[9741.218141981712, 13331.986680084825]\n"
     ]
    }
   ],
   "source": [
    "print(makeConfidenceInterval_std(data, 80, 0.95))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 置信区间验证\n",
    "由于置信区间不一定会包含要求的值，因此接下来，对于置信空间的准确性进行检验验证：     \n",
    "首先取出数据集的实际方差与标准差，然后多次计算置信区间，记录置信区间包含实际值的频率。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.843\n"
     ]
    }
   ],
   "source": [
    "pop_var=data.var()\n",
    "times_in_interval=0;\n",
    "times = 10000\n",
    "for i in range(times):\n",
    "    interval=makeConfidenceInterval_var(data, 80, 0.95)\n",
    "    if (pop_var>=interval[0]) & (pop_var<=interval[1]):\n",
    "        times_in_interval += 1\n",
    "        \n",
    "print(times_in_interval/times)   #置信区间有95的可能性包含总体方差"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.836\n"
     ]
    }
   ],
   "source": [
    "pop_std=data.std()\n",
    "times_in_interval=0;\n",
    "times = 10000\n",
    "for i in range(times):\n",
    "    interval=makeConfidenceInterval_std(data, 80, 0.95)\n",
    "    if (pop_std>=interval[0]) & (pop_std<=interval[1]):\n",
    "        times_in_interval += 1\n",
    "        \n",
    "print(times_in_interval/times)   #置信区间有95的可能性包含总体标准差"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "遗憾的是，似乎自己写的算法还存在一些问题，实现不了预定的置信度，如果直接把均值估计的代码拿过来呢？     \n",
    "下面使用的是均值估计的代码，唯一的不同只有将统计数据从均值换成了标准差："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "def makeConfidenceInterval_std_new(data, sample_size, _alpha):\n",
    "    sample=np.random.choice(a=data,size=sample_size)\n",
    "    sigma=sample.std()/(sample_size**0.5)\n",
    "    return stats.t.interval(alpha=_alpha,    #置信水平confidence level\n",
    "                 df=sample_size-1,  #自由度Degrees of freedom\n",
    "                loc=sample.std(),\n",
    "                scale=sigma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(8083.782109409111, 12711.55845781877)\n"
     ]
    }
   ],
   "source": [
    "print(makeConfidenceInterval_std_new(data, 80, 0.95))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9313\n"
     ]
    }
   ],
   "source": [
    "pop_std=data.std()\n",
    "times_in_interval=0;\n",
    "times = 10000\n",
    "for i in range(times):\n",
    "    interval=makeConfidenceInterval_std_new(data, 80, 0.95)\n",
    "    if (pop_std>=interval[0]) & (pop_std<=interval[1]):\n",
    "        times_in_interval += 1\n",
    "        \n",
    "print(times_in_interval/times)   #置信区间有95的可能性包含总体标准差"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "效果有所改善，不过感觉这并不是一个正确的方法，本来不准备写在作业里，但是表现不错实在舍不得删掉。   \n",
    "也许是因为本题中均值与标准差数量级相当，因此才有较高的准确性？      \n",
    "P.s. 本来想调用MATLAB计算的，但是树莓派装MATLAB好像很麻烦，所以最终没有。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 假设检验\n",
    "由于之前的数据生成环节是随机截取一个年份的数据，因此在这里就不太好假设均值数据，因此采用这样的策略：     \n",
    "进行两次不同的假设，第一次的假设值是实际平均值加上50，模拟猜中了的情况，第二次的假设值是实际平均值加上5000，模拟猜错了的情况。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. 假设正确\n",
    "首先是模拟假设正确的情况，"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "总体的均值为:  8421.026851851851\n",
      "空假设H0是：总体的均值是 8471.026851851851  \n",
      "\n",
      "从样本构造的t统计量 =  1.6940948692901128\n",
      "单样本双边检验的 p =  0.0941865267034535\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xd8VfX9+PHX+44MEhIgEFbCCFO2TBkVZ0W0qNUqrbVqHXWg9Wete9daqtaqX1eVWqvVInVSq6iIAwVFluyRhBUCJIQMErLuPZ/fH0kwhIybcO5+Px8PJffec895J59z3vdz3+dzPkeMMSillIosjmAHoJRSyn6a3JVSKgJpcldKqQikyV0ppSKQJnellIpAmtyVUioCaXJXSqkIpMldKaUikCZ3pZSKQK5gbbhzcrLp061bsDav/Ki0rAjTqUOww4hqcqCIxARtg0i0YsuW/caYLi0tF7Tk3qdbN5b/7W/B2rzyoyXL3qbqwp8GO4yoFjPvbSaN1zaIRHLyyTt8WU7LMkopFYE0uauI8dXHX7F963a/rf/2y2+ntKTUb+tXyk6a3FXE+Orjr9ixteVvrF6Pt03rn/2P2SQmJbbpvUoFWtBq7iqy3X313eTtyaOqsorzLzufn/ziJwAs+2IZcx6dg9frJblTMo+/9jjlZeU8df9TbF67GUH41W9/xdQzp/Ldl9/x8hMvU11VTY9ePbjt0duIT4hn5pSZnHTWSSz7YhmxsbHc/eTdFBYUsmThEr7/9nteffpVHnjuAXr27nk4ntm3zCYmNoat67cybOwwfn3zr3nq/qfYtnkbXo+XS397KVN+PIWK8gpm3zKb7Vu2k56Rzv59+7npwZsYNGIQM6fM5G/z/0Zyp2TmzZnHh//5EICzLjqLC359AXtz9nLbZbcxfOxw1q1cR5euXXjoxYeIjYsNShu0RbUIOQkJVDidwQ4l6sV5vaSVleFu47TsmtyVX9z6yK0kdUiisqKSa865hhPPPBFjGR674zGefONJuqd3p6SoBIBX/u8VEton8NKClwA4WHyQ4gPF/Ovpf/HYvx4jvl08/37+38z7+zwuvfFSABLbJ/LSgpf46K2PePoPT/Onv/+JSadNYuIpE5k6fWqjMeXvzefpt57G6XTy4qMvMnriaG575DZKS0q59pxrGTNlDPP/NZ/2ye15+ZOX2bZ5G1eedeVR69m8djML3lzAs+88CwauPe9aRk4YSfvk9uRsz+Gep+7hltm3cP/19/Plh19y+nmn++mvbL+chATad+1Kn+RkRCTY4UQtYwwFxcXk7NtH39K2lQJ9Su4iMg14EnACc4wxs5tY7gLgP8A4Y8zyNkWkIsLbL7/N4o8WA5C/J5/d23dTVFDEiPEj6J7eHYCkDkkArPh6Bfc+de/h97ZPbs/ST5eyPXM7N1xwAwCeag9DRg85vMwpM04B4NQZp/LsQ8/6FNPU6VNx1vZIly9ezpKFS3jjxTcAqKqqIi83j7XL13L+5ecD0HdQX/oN7nfUetYuX8uUH08hvl08ACeecSJrlq1h8umT6Z7enf5D+gMwaPgg9ubs9Sm2UFHhdGpiDwEiQkpyMvn797d5HS0mdxFxAs8ApwM5wHciMt8Ys6HBcu2BG4Fv2xyNiggrs3JYsXIrz7z9DHHxcdw08yaqKqswxiA0kjQMRyUTYwxjp4zlnqfuaXQbRyzvYx6Kj48/Yv0PPPsAvfr1Omq7LWpmEXeM+/DPDocDr7dt9f1g0sQeGo61HXw5oToeyDTGZBtjqoC5wDmNLPcH4BGg4pgiUmGvrKKK9kntiYuPY2fWTjasqukHDB09lO+Xfc+eXXsADpdlxv5oLO+88s7h9x8sPsiQ44ewbsU6dm/fDUBFeQW7sncdXuaz9z87/O/Q44cC0C6hHYfKDvkU47gTx/HOP985nMy3rt8KwPCxw/n8f58DsH3rdrI3Zx/13hHjR/D1x19TUV5B+aFyFn+8mBHjR/j2x1EqQHxJ7j2BXfUe59Q+d5iIHA+kG2Peb25FInK1iCwXkeX5xcWtDlaFhwmDeuP1ern0tEt54c8vMOT4mnJKh5QO/O7h33HvNfdyxZlX8OCsBwG4ZNYlHCw+yOVnXM4VZ17BqqWr6JDSgdsevY0//PYPXDHtCq4/73p2Zu08vI2DxQe5YtoVvPXyW1x/z/UAnPyTk3njhTe46qyr2L1jd7Mx/uqGX+HxeLjizCu4/IzLeenxmnr/OZecQ1FBEZedfhkv/eUl+gzoQ0L7hCPeO3DYQM644AyuPfdarjvvOs666CwGDB1g29/PF6Nm3sSomTcFdJuBUFRUxLMvvujTsvc//DCPPfUUAPc+9BALP/usyWXfff99Nmza1OTrz//977zy+usAnDR9OstXrmxzzLl79nDBJZf4/H5/kZa+horIz4AzjDFX1j6+BBhvjLmh9rEDWARcZozZLiKfA7e0VHMfO2iQ0StUI8/CwhWs37SYkT/3X+KpP2rFLvklFq8srqK03IPl9eJ0x+Ao3cOav93Ja5+9ekS5JRTUJfbVc5846rUVuStwfbGY305sfRtsTE7muP79jzm+ttq+YwdnX3gh675tubp7/8MPk5iYyC033tjispddcw1nT5vGBeeee9RrHo8Hl+uHCvVJ06fz2EMPMXb0aNtjbq2NmZkc16AjLCefvMIYM7al9/rSc88B0us9TgNy6z1uDwwDPheR7cAJwHwRaXHjKvI8tONV/ukJr3PpldWGpz6qJHOflzip4vtnfs/yx69n8TN/oMdZ14EzvAaVvbom/Nqgzu333UfWtm2MmjyZ399991Gv//HRRxl4/PFM+fGP2bx16+HnL7vmGt58993D6xgybhwjJk7klrvuYsm33zL/gw/4/T33MGryZLKyszlp+nRuuu02xk6dypPPPXfEtwCAV+fOZdTkyQybMIFly2v+lg2XGTZhAtt37Dgq5u07djBswgQAKioquPzaaxl+wgkcP2UKn335JQAvv/YaP734Yqaddx4DRo3i1nsaP7d0LHzZa78DBohIX2A3MBP4Rd2LxphioHPdY1977kq11dyv5tq2LmMMf/+8kpwCi/83PZYRvRLg/BcAWJbl4dlPKvnXV1VcdmJM1J1ofODjLDbss/eK3CFdE7nvx0ePQKoz+4EHWLdxI6u//vqo11asWsXct95i9ddf4/F4GP2jHzHm+OOPWObAgQO889//smnFCkSEoqIiOnTowIzp04/quVdVVbH8iy+AmsRd36HyclZ//TVffv01v77++mZ75Q1j3r7jhwvpnqkt16z95hs2bdnCj889ly21JZ/Va9eyavFiYmNjGTRmDDf85jekp6U1uZ3WarHnbozxALOAj4CNwDxjzHoReVBEZtgWiVJB8OHqapZleblggpsRvY7s64zv5+Ls4918sdHDZxs8QYpQ1Vm8ZAnnnX027dq1IykpiRnTpx+1TFJSEnFxcVw5axZvz59Pu3btmlzfReef3+RrP7/gAgBOnDyZkoMHKSoqalPMXy1dyiUzZwIweOBAeqensyUzE4BTp04lOTmZuLg4hgwaxI5du5pbVav59H3TGPMB8EGD5+5tYtmTjj0spfxvzU4P//m2mvH9nEwf1XhN/afj3OwqsHjt6yp6dnIwqHv0XLnZXA87WFr69uRyuVj22Wd8+vnnvPneezz9wgsser/xcR4JzST+htsREVwuF5ZlHX6uoqLlgYHNndOMjYk5/LPT6cTjsbcDoXPLqKhUVml4fmElaSkOrjgptsmk4XAIvzk1li7thWc+rqSyum2XgivftE9M5GATV2SeOHky77z/PuXl5Rw8eJD/fvjhUcuUlpZSXFLC9DPO4K9/+hPfr13b4nob88bbbwM1Pe/kpCSSk5Pp06sXK1evBmDl6tVsqy2/tBTza/PmAbBl61Z25uQwaEBgRlaF15kiFfL+NvBmVq75JNhhtOirTR4OVcGtJ8UQ626+N9guVrhsaiyz51fwTaaHqccFf+TM5j/e3ORrN0+8GfeHod8GjUlJSWHyhAkMmzCBM08/nUcfeujwa6NHjeKin/6UkZMmkdqlC+MaGc1ysLSUc2bOpKKyEmMMj9fW0mdecAFX3XADTz3/PG++8kqLccTFxnL8lClUV1fz0jPPAHD+Oefwyr//zdDx45kwdiwDa0cVNYz5+quuOrye6668kmtuuonhJ5yAy+Xi5eeeIzY2MHMNtTgU0l90KGTkCvWbdVjGcMe/y2kfL9x9XnzLb6Dm6/U9/ynH6RDuPz8u5E+utvVmHcEeCqmO5O+hkEr57L/7l/CVd1uww2jW+hwv+0oMpw7zvQcuIpw6zM2O/RZZ+6yW3+BnKQuXkLJwSaOvLdkV+m2g/E+Tu7LVX3LmMdezKthhNGvhOg9J8TA2o3UnRycOcBEfAwvXVfspMt+lz5lH+px5jb42b33ot4HyP625q6iSX2KxZoeXs0e7cTtbV1qJcwtTBrpYtMHDzw8ZktsFrzSz/tkHgrZtFR60566iymcbPIjAyUPa1q85ZZgbrwVfbAxu7726UzLVNk6/oCKPJncVNao8hi82VjO6r5NOiW3b9bt3cDA0zcHnGzx4reANi+z25gK6vbkgaNtXoU+Tu4oa32Z6KKuEU4ce21DGU4e6OVBmWLU9eHO1a3JXLdHkrmz16uA7uccdmreVW7TeQ4+OwuAex7bbj+rtJCVRWLQ++CdWG3Pnj0K3Dez28muvMet3v2vTe+tPNtbc+nP37GnVeutPHNbw+dfnNX4S3B80uStbpcel0tXRPthhHKXgoMW2fIspg1zHfocbhzBpoIuNuRalFaF3xWpqQmi2QThqS3JvyvadO3n9P/+xZV2+0OSubPVG3iIWerYEO4yjfL+zpoQyqrc9A8SO7+3EGFizM/Ruo7doW2i2ga9eef11RkycyMhJk7ik9mrP/374IRNOPpnjp0zhtBkz2JeXd9T79uXlcd4vfsHISZMYOWkSS7799qhe9GNPPXXUDJAAD86ezbipUxk2YQJX33gjxhjefPddlq9axcVXXsmoyZMpLy9nxapVTD3zTMaceCJnnHsue/bW3CN3xapVh7f7TBM3G7n9vvtYvHQpoyZP5q9PP23Hn6pZOhRS2eq53PmUePdzYrADaWD1di9dk4TuHewZvtgn1UFSvLB6h4dJA0PrMJq/eT4O734andmvtRqZefEo06ZB3Q0zpk+Hiy+u+a+gABrekeiDD45+fz3rN27kj489xteffELnlBQOHDgAwJQTTuCbRYsQEeb885888sQT/KVBkr7x1luZOnky77z+Ol6vl9LSUgp9nM1x1tVXc+/ttwNwyVVX8f6CBVxw7rk8/cILh2/cUV1dzQ2//z3vzZ1Ll86deeOtt7jrwQd56dlnufy66/i/Rx9l6pQpjc5DDzVTAz/21FO8H6Dee2jtlUr5QWW1YUOul1OGHHtJpo5DhJG9nSzP9uDxGlytHDOvGrfoiy+44Nxz6ZySAkCnTp0AyMnN5aLLLmPPvn1UVVXRt3fvRt/7Su2UJk6nk+TkZJ+T+2eLF/PIE09wqLycA4WFDD3uOH5y5plHLLN561bWbdzI6efU3ELa6/XSvWtXiouLKSouZuqUKQBcMnMmH34S/Ll9NLmriLc+x4vHC6P62Lu7j+rtZPEmD1v3WhzXM0KnAm6hp93s8ikprX6/MYbGPiZv+P3vuXnWLGZMn87nixdz/5/+5NP6fJmmt6KigutuvpnlX3xBeloa9z/8cKPLGWMYOngwSz/99Ijni4qKQnKuIa25q4i3eoeX+BgY2M3e3X1omhOXA1bv0Bt52OXUk05i3jvvUFBQAHC4LFNcUkLP7t0B+GftjayPeu/UqTw3Zw5Q06suKSmha2oqefn5FBQUUFlZyfsLjh4+WpfIO6ekUFpaypvvvXf4tfrT+Q4aMID8/ftZWntXpurqatZv3EiHDh1ITkriq6VLAQ5P8dtQa6cdPlaa3FVEs4xhzU4vw9KctpdO4tzC4J5Ovt8ReidVw9XQ447jrltuYer06YycNImb77wTgPvvuIOfXXopY0488XDJpqEnH3mEzxYvZvgJJzDmxBNZv3Ejbrebe2+7jfEnn8xpM2YweODAo97XoUMHrrr0UoZNmMAZ5513xFTCl118MdfcdBOjJk/G6/Xy5quvctt99zFy0iRGTZ7MktpE/49nn+X63/2OUZMnN3mDjhHDhuFyuRg5aVJATqjqlL/KVvuri1m28n3aXXRxsEMBYFuelwferuCqU2KYPND+edgXrqvmX19VMXtmPN06BK6v5D5QMw1sY1MQFFcU437vfaZPaH0b6JS/oUWn/FUho7M7mQ7i2xzpgbB6hxcBRqT75/TSyF7Ow9sJpObmlkmOC602UMGhyV3Z6uW9C/ifZ2Owwzhs9Q4v/bo6aB/vnxNeXZIcpHWSgNfdm5t+YEFmaLWBCg5N7spWL+9dwIfe0EgshaUWO/ZbjOrj35Eso3q72LLHoqwycCXOlpL7sbRBsEq16kjH2g6a3FXEsvuq1KaM7O3EMrB2V+BKM6vnPsHquU/Yvt44r5eC4mJN8EFmjKGguJg4b9v3KR3nriLW9zu9pCQKPTv6dwxyv1QH7ePg+x0eTugf3odUWlkZOfv2kb9/f7BDiXpxXi9pZWVtfn9474lKNcGyDJtyvYzvZ99VqU1xOITjejrZuNuquQgnABe0pL/wBgC7rr7I1vW6jaFvAMdiK//RsoyKSLsOWJRXwaDugblydHAPJ0WHDPklgSlnpCxaSsqipQHZlgpP2nNXtvpg+Gy+Wf5eywv62abcmkvOBx3j3O2+qvsQ2bTHS2pycPtMs0+bTcxbwW8DFVzac1e2aueMI07sv1iotTbneunSXkhp4+30WqtHR6F9HGzOtVpe2M/iXKHRBiq4NLkrWz27+13e9qwJagyWMWzZ62VQj8BN5iUiDOzuZPOe4E9F8O6m4LeBCj5N7spW8/I/Z5E3M6gx5BYaSitgUPfA7t6DezjZf9BQcDC4vffPtwe/DVTwaXJXEWdTbk3veXAAe+7ww4fJphDovSulyV1FnC17vHRKEDq3D+wc22mdHLSLCY26u1Ka3FVEMcawKddiUA9HwG+g4HCETt1dKU3uKqLsLTaUlJuAjW9vaFAPJ/uKDUVl2ntXwaXj3JWtPh/1BEuWvU1VkLa/ubbeHsiRMvX9UHe3OKG///pOzc0r88S0J4iZ97bftq3Cg/bcVUTZvMdLUrzQLTk497Ts3dlBnBu25GppRgWXJndlq8d2vcHr1SuDsm1jDJtzLQYHod5ex+kQBnRz+n3ETPoLbxyeX6ahN9YFrw1U6NCyjLLV+wVLKbH2c0EQtp1/0HCgLHj19jqDejh481svJeWGJD/dJCRp1fomX1uasxSHpbM6RjtN7ipiBLveXqfmw6WaLXu8jM3wzyG2/rkH/bJeFTm0LKMixpY9FolxNfO8BFPfLg5iXOiQSBVUPiV3EZkmIptFJFNEbm/k9WtEZK2IrBaRr0RkiP2hKtW8rDwv/bo6cQSp3l7H5RT6dnGQtc9/wyH7PvIifR950W/rV+GvxeQuIk7gGeBMYAjw80aS9+vGmOHGmFHAI8DjtkeqwkK8I5ZYAl8WKas05BYa+qWGxpfRfl2d7NhvUeXxz/zuySvXk7yy8bp7rCs4baBCiy9Hwngg0xiTbYypAuYC59RfwBhTUu9hAqA3YIxSH474M3+JPaflBW2WnVdTAunfLTSSWr+uDrwW7Nwf+IuZ/nxacNpAhRZfkntPYFe9xzm1zx1BRK4XkSxqeu43NrYiEblaRJaLyPL84uK2xKtUo7L2WQg19e5QUPcNwp+lGaWa48uR0FgB86ieuTHmGWNMP+A24O7GVmSMecEYM9YYM7ZLcnLrIlVh4Q87XuEf1csCvt3MfRY9OwnxMcGtt9fpkOAgJVHI3Bf4k6qvfB+cNlChxZfkngOk13ucBuQ2s/xc4NxjCUqFr08LV7LCygnoNi1jyN5XczI1lPTv5iArL/A995V7At8GKvT4kty/AwaISF8RiQFmAvPrLyAiA+o9PAvYal+ISjVvb5HhUBX07xoaJZk6/VKdHCg1FJZqaUYFXotXWBhjPCIyC/gIcAIvGWPWi8iDwHJjzHxgloicBlQDhcCl/gxaqfqyaksfodZz71f7YZOZZzEuQPdyVaqOT5fPGWM+AD5o8Ny99X7+rc1xKeWzrH0W8THQrUNo1Nvr9O7swOWsiW9cRrCjUdFGpx9QtkpxJwGlAd1mVp5Fv9TgX7zUkMsp9OnsINsPJ1WrOyY1+VpSbBKOALeBCj2a3JWt3hr6YEDncy+vMuQcsBg92h2gLbZOv64OFq334PEaXE77Pnyam1vmwZMf1Pnclc4to8LbtnwLY36ob4eafl2dVHthV4GeVFWBFZpHhApbd2S/yHPVSwK2vbqTqRmpoXUytU7dh47dQyKbm1vmxRWBbQMVmjS5K1stLVnPemtvwLaXtc+iWwchMS606u11OiUIHdrJ4Q8hu7gLi3EXNn6V9/r8wLaBCk1ac1dhyxhD1j4vI3uH7m4sIvTr6iDT5mkItvzpFlvXpyKP9txV2MovMRysCN16e53+XZ3klxhKynU+PRU4oX1UKNWMut5wqEzz25TDdXcbSzMD73iMgXc8Ztv6VOQJ7aNChZ202C50kcSAbCs7z0uMC3p2Cu3duE8XBw6BbBtPqrbblkO7bY3PH9MlIXBtoEJX6BYrVVj613F3BWyce3aeRd8uDpyO0DyZWifGJaSnOA7POe9vd/3oLmL26Dj3aBfaXR6lmuDxGnbut+gbokMgG+qb6mBbnoVltO6uAkOTu7LVTZlP80TVl37fzq4CC48FGSFeb6+TkergUBXkFfs/uT+9LDBtoEKblmWUrVaXZlJi9vt9O3X16/BJ7jXfMLLzLLp18G/MmQcycQSgDVRoC48jQ6kGsvMskuIhJTG06+11enQQYl0ErO6ulCZ3FZa25Xvpm+pEQmwmyKY4HEKfLg5bR8wo1RxN7irsHKo07Ck0YVOSqZOR6mTnfguPV0+qKv/Tmruy1cB2aewr9fh1G9v3WxjCp95eJyPVgceqORl8rKN8DvVNa/K1tKQ0nPv92wYq9GlyV7Z6YeAtLCny7zj3urp13y7hMQyyTt2HUXbesSf35uaWuWXSLcTk6Dj3aBdeXR+lgG15Fl2TQncmyKZ0ShSS4kXr7iogNLkrW1295TH+XLXIr9uo6fmG364rImSkOtiWf+wjZpqbW+axJf5vAxX6tCyjbLXlUA4lpshv6y8stSgsMyF7c46WZKQ6+H6Hl0OVhnaxbf/mUd0xucnXckpycPixDVR40OSuwkp2fnhdvNRQRqoDQ81J4SE92/4Bte3Wq+wLSkWk8DxCVNTalmfhdECvzuG56/bpUnelql7MpPwrPI8QFbWy87ykdXIQ4wqvk6l1EuOErsnCtmM8qTr02nsZeu29NkWlIpGWZZStRiX2J7fUP+u2jGFbvsUJ/cN7t81IdbAp99iSu7uwpMnX+nfqj0Onlol62nNXtnqi/yxuijnRL+veW2QorwrfenudjFQnhWWGwlL/DImcNd5/baDCR3gfJSqqbKu7eClMR8rUqRvGWXdyWCl/0OSubPXLjX/kgaqP/bLurDyLOHfNDIvhrFeKA6fD3tvu1ffHxf5rAxU+wrt4qUJOTmU+JcY/Rfe62+o5Qvy2ei2JcQm9Uhxk23jD7Pryy/Jx+KkNVPjQnrsKC1Uew64Ci4yu4V2SqVNzpaqFZekMkco/NLmrsLBzv4U3jG6r15KMVAcV1bCnSJO78o/IOFJUxMsKs9vqtaTuG0iWXsyk/ERr7spWE5OGklO22fb1bsvz0ilB6JgQGcm9a7LQLqbmPMKJg1v//uLRQ5t8bWiXoTgK7G8DFV40uStb/SnjKpbst38+96w8i4yukZHYARwi9E11kL2vbSNmmptb5qoxVxGTpfO5R7vIOVpUxCopN+SXGDK6RNbumpHqJOeARWW11t2V/SLraFFBd/76e7mz8gNb11l38VKkjJSpk5HqwDKwY3/re+/NzS1z72f2t4EKP1qWUbYqqC6hhApb15mdZyECfSKs597v8ElVi4HdW/fBVXJ80zX3ksoSHDa3gQo/mtxVyMvKs+jZUYhzh/fFSw0lxQud20vtxUzuVr1319UX+ScoFTEiqyukIo4xhm153rC981JLMlIdek9V5Rc+JXcRmSYim0UkU0Rub+T1m0Vkg4isEZFPRaS3/aGqaLSv2FBWCf0iaKRMfRmpTgpKDcWHWndSddTMmxg18yY/RaUiQYtHjIg4gWeAM4EhwM9FZEiDxVYBY40xI4A3gUfsDlSFh1M7jmaMI8229WUfvngpMnvu/epmiLTxYqbR3e1tAxWefOkOjQcyjTHZxpgqYC5wTv0FjDGfGWMO1T78BtA9K0rd0/tXXO4eb9v6svO8xLqgZ8fIqrfX6dXZgUPsnSHyVyPtbQMVnnxJ7j2BXfUe59Q+15QrgA+PJSil6mTnWfSJgJkgmxLrFtL9OEOkil6+JPfGjqpGC4Qi8ktgLPBoE69fLSLLRWR5fnGx71GqsHHmmtv4XeV7tqyr2mvYud+K2JJMnYxUB9n5Fpax52Km2xba1wYqfPmS3HOA9HqP04DchguJyGnAXcAMY0xlYysyxrxgjBlrjBnbJTm5LfGqEFduVVKJPb3QnfstPBYRNe1AYzJSHZRXwZ5Ce5J7pce+NlDhy5ej5jtggIj0FZEYYCYwv/4CInI88DdqEnue/WGqaJS5t6YO3T/Ck3v/bjXfTDK1NKNs1OJRY4zxALOAj4CNwDxjzHoReVBEZtQu9iiQCPxHRFaLyPwmVqeUz7bu89K5feTMBNmUbslCQixktnESMaUa49MVqsaYD4APGjx3b72fT7M5LhXljDFk7rUY3COyEzuAiNC/m5PMvdpzV/bR6QeUrc5Omcj2srXHvJ6CUkPRIXO4ZBHpBnR18P0OL6UVhsS4lkcGFZwyscnXJqZNxFlw7G2gwpsmd2WrW9IvYske9zHP515XohjQLfJ77lBXd68ma5+Xkb1bPiybm1vmomEXEbOhdXPVqMiG2ZfBAAAWlklEQVQTHUeOCjuZe2suXkrrFB27aN8uNRczad1d2SU6jhwVMCetvolZlcd+F6DMfRYZqQ6cEXrxUkOxbqFXZ4fPdffm5pa5aYE9baDCm5ZlVMiprK65eOms46OrtNC/q4MvN3nwWqbFD7W9F0wLUFQqXGlyVyEnO8/CMpE/vr2h/t2cLFznYVeBRZ8uzZ9I1uSuWhJdR48KC3UX8/SLsNvqtWRA7YeZL3V394Fi3Ad0Cg/VNE3uKuRk7rXo0VF8GhIYSTolCh0TxKe6+9Dr7mPodfcFICoVrrQso2x1YZeTyD60us3vt4whc5+XMX2jb9cUEfp3dbB177GNmDmpz0m4DrS9DVRkiL4jSPnVdT3PZcluq83j3PcW1dx5qX+UjG9vqH83J99leyksteiY2La/wbmDzyVmjQ6pjHbReQQpvznkraDCVLf5/XX19v5RVm+v078VdfemVHiOrQ1UZNDkrmw1fe3t3FL13za/P3OvRUIsdOsQXfX2Or07O3A5j22GyNsXHlsbqMigyV2FlMx9Xvp1deKQ6EzuLqeQ0cWhV6qqY6bJXYWM0gpDbqGJuvHtDfXv5mR7vkVltT0371DRKbqPIhVSNu+pKUUM7hGd9fY6g7o78FqQpb13dQw0uauQsSnXS4wL+qZG9245sLsTEdi0R+d3V22nQyGVrS7rNo2th1a06b2bci36d3XgdkZnvb1OfIzQp7ODTblNJ/fmph+Y1n8azgNtawMVOTS5K1td1m0aS3YeavU499IKQ06BxXnjomuysKYM7uHkk7XVVFYbYt1Hf9i1lNxjVh7yZ3gqDET3919lu/3VxRSZ8la/b/MeLwatt9cZ3MOBx4KsvMbr7s3NLVNc0bY2UJFFk7uy1QXr7+Puqg9b/b5NuV7cTq231zlcd2+iNNPc3DL3fd62NlCRRcsyKiRszrXo303r7XXq6u6bm0juu668MMARqXCjyV0FXWmFYVeBxblabz/CoB5OFq6tpspjiHEd+aFXcNqkIEWlwoV+B1ZBt0Xr7Y06XHdvZLx7fNZO4rN2BiEqFS40uaugq6u3Z2i9/QgDuzVddx901+MMuuvxIESlwoWWZZStru0xg82Zy1r1ns17dHx7Y9rFCr1bGO/emBmDZuA60Lo2UJFHu0rKVhelnsJproE+L19WWXMzbC3JNG5wDwdZeRZVHt/nmTmlb+vaQEUmTe7KVrsq8thnHfR5ea23N29wDyceb+vmmckra10bqMikyV3Z6pJND/OH6k98Xl7HtzdvQDN196Y8vLh1baAikx5RKqjq5pNpONRP1UhoY91dKU3uKmiKDxl27Lc4rqeWZJozpKeTzH0W5VU6v7vynSZ3FTTrdnkAGNFLk3tzhqc78VqwYbf23pXvNLmroFmzy0tSvNCrs+6GzRnQzUGcG9bu1OSufKfj3JWtfpd2IRu3Lm1xOcsyrNvlZVRvV9TeL9VXLqcwNM3Jml1ejDGISLNzy1w49EJcX7XcBiqyaXJXtvpJ50mkZO9tcT737DyLskoYriUZnwxPd7Jim5fcQkPPTtLs3DKT0icR49wbwOhUKNLvw8pWmw/tZIdV2OJya3d5EYFhaZrcfVH3IbhmV01pprm5ZXYW+9YGKrJpcle2+s2Wx3m0+rMWl1uz00u/VAeJcVqS8UVKooOeHYW1O2tOQjc3t8zjS31rAxXZtCyjAq6k3LA9X6f4ba3hvVwsXFtNRbUh+/dXBjscFeK0564Cbt2umikHRqRrSaY1RvRy4rFg424vJWOGUTJmWLBDUiFMk7sKuDU7PSTFQ+8uuvu1xsDaIZFrdnpJWrGOpBXrgh2SCmE+HV0iMk1ENotIpojc3sjrJ4rIShHxiMgF9oepIoVlGdbleBmWrkMgW8vlFI7r6WTNTi99H51DxqNzgh2SCmEtJncRcQLPAGcCQ4Cfi8iQBovtBC4DXrc7QBVe7u59CZe6xjb5+rZ8i9IKLcm01YheTgpKDVWeppe5ZETzbaCigy8nVMcDmcaYbAARmQucA2yoW8AYs732Nd/nJVUR6bSOY2jn3NHkOPc1O70IMFSHQLbJ8NoPxbIKQ2xi4998xvQYQ4xzRyDDUiHIl7JMT2BXvcc5tc8pdZTVpZlssfKbfn2Hl4yuDtrHa0mmLTq3d9Cjo1Ba2fQkYpkHmm8DFR18Se6NHYVtmp5ORK4WkeUisjy/uLgtq1Ah7qbMp3mqenGjr+UVW+zYbzE2Q0fgHotxGS4OVRo8TXxPfnpZ022goocvyT0HSK/3OA3IbcvGjDEvGGPGGmPGdklObssqVBhbllVTKB6XoSWZYzG+X82H48FynQJYNc2X5P4dMEBE+opIDDATmO/fsFQkWpblpV9XB53b6xDIY9GzU83NTUo0uatmtHiUGWM8wCzgI2AjMM8Ys15EHhSRGQAiMk5EcoCfAX8TkfX+DFqFn71FFjsLLMZrScYWSfFCeZWhsEzHMKjG+XSkGWM+AD5o8Ny99X7+jppyjVKNOlyS6aclGTu0jxf2H4Tl2V5OH67fhNTRtBulbPVw3ytZu+GLo57/LtvLgG4OOiVqIrLD7juu5OVPK8jO8nD68CPn6Lly9JW4Pz26DVR00eSubDUpeRg4txwxzj230GJXgcXFk2OCFlekKRkzjHZUseW7ag6UWkd8aA5LHUaMc0sQo1OhQLtRylZLitex1rvniOe+y/IgwNi+WpKxS9KKdZx5cDNQ862ovnV5R7eBij6a3JWt7tw2h795jrzF27IsDwO6O+ioJRnbZDw6h3HP/530FAffZR05F8GclUe3gYo+WpZRfrX7gMXuQsMvp+jc7Xba/MebAZhQ7OTNZdUUHLRI0SGmqh7dG5RffVtXktELl2xV3q8X5f16Ma72gqZvs5qZSUxFJU3uym+8lmHxJg9D0510aKe7mp1SFi4hZeESuiY76N/VwZcbPRijFzWpH+gRp/xm9Q4vhWWGU4Zo9c9u6XPmkT5nHgAnD3Wxt9iwcbde0KR+oMld2eqJ/rO40f0jABatr6ZTgjCyt5Zk/GlchovEOFi0oRqAWeN/aAMVvTS5K1uNSuzPQEcX9hZZrM+xmDrEhdOh0/v6U4xLmDLIzcptXgrLLPp3qmkDFd00uStbLSxcwXfenXy+sRqnA6YO1pJMIJw8xIVl4MuNHlbk1rSBim565ClbPbTjVYqq91O9ycPxfZx0SND+QyB0TXYwLM3J5xs9WPtfxenZz2+DHZQKKj3ylO1KTAxllXDqUB3bHkinDHVRWGYordBRM0qTu/KDQiuW7h2EwT109wqkkb2ddEoQiso0uStN7spmZV5DuXFx8hA3InoiNZCcDmHqEBdllYZKo4d2tNM9QNkqt9zgwDB5kJ7OCYapg12IQIEVF+xQVJDpEahss/GgBWXX87PYnSTEaq/dn+rmlmmoQ4KDs3v/lqVbq9leZtFHT2hHLW15ZZvHt1aR4kzjl3HlwQ4l4tXNLdOYi0/IoJ3pyZNZVY2+rqKDJndli9VFXj7J8zKx2wpWma3BDifi1c0t05gNBd8wMO5/vJvrZWupTkkQrTS5K1v8ZWs1ndywofpt5npWBTuciFd/bpmG5q2fx07HByQ44a9btfcerbTmro7Ztwe8LC7wctegGF4rDnY00WH9sw80+7pLDL/s4+aprGrWlXgZlqTz+0Qb7bmrY2KM4S9bq0iNFX7ZS/sKgVLdKZnqTsnNLnNFHzfJbvjr1uoARaVCiSZ3dUwW5ntZVmgxK8NNvFNHyARKtzcX0O3NBc0uk+wWru7j5tN8L0sKvM0uqyKPJnfVZsXVhrvWVzE4UZiZrr32QPIluQP8uo+bPu2E29ZVcsijV65GE03uqs0e2FhFQZXhseGxxNRO6/vq4Du5x316kCOLbnf+6Ic2iHcKjwyLJafc8OctenI1mmhyV22yMM/D27kers9wMyz5h5N16XGpdHW0D2JkKjXhyDYY38nJZb1d/HOnh6VanokamtxVqxVXG+6sLcfM6nfkzI9v5C1ioWdLkCJTAIu2Hd0Gtw6MoXc74VYtz0QNTe6q1Rorx9R5Lnc+73rXBSkyBTB/89FtEO8UHtXyTFTR5K5a5fVd1byd6+G6BuUYFfrql2fey/UEOxzlZ5rclc8+zfNw9/oqTurs5MZ+eiOOcHTbwBjGd3Rwy9pKHR4Z4TS5K5+sLvIy6/tKhiY5eGZULG696XVYinMKL46Oo2+C8JtVFWw6qHPPRCpN7qpF28ssrlhZQecY4aUxcSS4NLGHs2S38I8xcbRzCpctr2BPhSb4SKRXnqhmZZdZXLa8AsvAP8fG0aWFedrfHPoAy1a+H6Dooldzc8s8cNIDuN9rvg16xjv4x5hYLvy2gku+q+AfY+JIb6d9vUiiramatHi/l3OXllPqMfxjTBwZPtz4obM7mQ4SH4Dooltzc8skx/nWBkOSnMwZE8e+SsO535TzXaHW4COJJnd1FGMM/9xRzWUrKugeJ7w3MZ5RHXwbGfPy3gX8z7PRzxGq5qYfWJDpexuc0MnJuyfEk+QSfrGsgnk5OslYpNDkro5QVGW4bV0V922s4uTOTt46Ib5VX9df3ruAD72a3P2tpeTemjbol+jg3YnxTOjk4NZ1VdyzoZKDeqFT2NOauwLAMoZ5OR4e2VJFUTVcn+Hm5gFunKInT0PR6rlP2Lq+upOsf9pcxT92eFiwz8udg2I4t7sT0X0gLGnPPcoZY/jmgJfzvqng9vVV9E908L/J8fx+YIwm9ijjdgj3HhfLuxPj6BEn/L81lVy4rIIVhV6M0Z58uNGee5Qq8xjeyfXwr53VbCo1dIkVnhgRyznaUwsL6S+8AcCuqy+yfd0jk528c0Ic/9nt4c+bqzj/2wqGJTn4VS8XP+nu0nn7w4RPyV1EpgFPAk5gjjFmdoPXY4FXgDFAAXCRMWa7vaGqY1VUZfiywMvn+V4+3ueh1AtDkxz8eZibGXrQhpWURUsB/yR3AIcIF6W5Oaubi3dyPby6s5pb11Xxx81V/DjVxUldnExJcZLs1n0mVLWY3EXECTwDnA7kAN+JyHxjzIZ6i10BFBpj+ovITODPgH/2OuUTj2XILjOsP2ixvsTLyiKL1UUWFtDBDT/u6uLidBejOzhs7al/MHw23yx/z7b1qdabfdpsYt6ypw0SXcIlvdz8Mt3FskKL13ZVs2Cfh//s9uAUGN3BwfHJToYmORia5KBvgmg5L0T40nMfD2QaY7IBRGQucA5QP7mfA9xf+/ObwNMiIkYLdbYyxnDIW1NSKfVCYZXhQLXhQJWhoMqwu9yQU27IKbfYVW6oqr3wMNYBQ9o7mNXPzUldnIxMdvjtAGznjCNO3Oi8g8ET54ojRuyd+0dEmNDJyYROTjyWYVWxxRf5NTdGf3lHNVW1R3qcA9LihbR4B+nthB5xQkqM0Kn2v45uIcEFiU4h3omWAP3Il+TeE9hV73EOMKGpZYwxHhEpBlKA/XYEWd/CPA/vBnBGO18+neo+wuqWrf+vMbX/1j5hAAuwTO1/GLy1P3sNVFtQbcBjDNUWVFpQ6TVUWnDI23w8Hd2QFu9gYKKDU7s4GNxeGJrkpF+C4ArQXDDP7n6XbM8azuanAdmeOtq7m97F5VnDJD+1gcshjOvoZFxHJ7cA1ZYhs8ywrtjLpoMWORU1nYxVxV6Kmxk27wDinTWdj1inEOsAtwNcIsQ4wCXgFHDU/YvgqH3sAISa/0ntz1L7MzTyrw+7fyA/Zi5Kc/Ojzv6dVdWX5N7Y79wwx/iyDCJyNXA1QK/OnWHvXh82f6SCAjcbCgM7I2Fzjd5w52ly56pdxkHtjlm7IzoFnBgcAm6gvYBLDG4HuF0Q5zDEOiBGoJ3DkOCABKch0WHo4DKkuAydXIZOTkO7xvaVstr/AmTe7o856KhgbGnr21b5rspb891obyN/54+zPibeUdGm46st3MBxwHEuoOORr5V54YBHOOAVDniEIo9QagllXqHUgnJLajowRqiwwGNqOjfVluAx4KWm41ON4DU/dJgsav6r6zxBM/+aljtpgS4xFMaWgse/nVRfknsOkF7vcRqQ28QyOSLiApKBAw1XZIx5AXgBYOzYsYaZM1sd8EVoMT+kvfw87YGZw1rftqoVEp4HGv87P7/8eUiAthxfdkuo/S+9pQWV737+c58W82Wc+3fAABHpKyIxwExgfoNl5gOX1v58AbBI6+1KKRU8Lfbca2vos4CPqBkK+ZIxZr2IPAgsN8bMB/4OvCoimdT02IPfZVBKqSjm0zh3Y8wHwAcNnru33s8VwM/sDU0ppVRbSbCqJ2PHjjXLly8PyraVUipcicgKY8zYlpbTuWWUUioCaXJXSqkIpMldKaUikCZ3pZSKQJrclVIqAmlyV0qpCKTJXSmlIpAmd6WUikCa3JVSKgIF7QpVEckHdti82s74YQ75ANPfIfjCPX7Q3yFU+ON36G2M6dLSQkFL7v4gIst9uSw3lOnvEHzhHj/o7xAqgvk7aFlGKaUikCZ3pZSKQJGW3F8IdgA20N8h+MI9ftDfIVQE7XeIqJq7UkqpGpHWc1dKKUWEJncRuUFENovIehF5JNjxtJWI3CIiRkQ6BzuW1hCRR0Vkk4isEZF3RKRDsGPylYhMq913MkXk9mDH01oiki4in4nIxtr9/7fBjqktRMQpIqtE5P1gx9IWItJBRN6sPQ42isjEQMcQccldRE4GzgFGGGOGAo8FOaQ2EZF04HRgZ7BjaYNPgGHGmBHAFuCOIMfjExFxAs8AZwJDgJ+LyJDgRtVqHuB3xpjjgBOA68PwdwD4LbAx2EEcgyeBBcaYwcBIgvC7RFxyB64FZhtjKgGMMXlBjqet/grcCoTdSRFjzMfGGE/tw2+AtGDG0wrjgUxjTLYxpgqYS01HIWwYY/YYY1bW/nyQmqTSM7hRtY6IpAFnAXOCHUtbiEgScCLwdwBjTJUxpijQcURich8I/EhEvhWRL0RkXLADai0RmQHsNsZ8H+xYbPBr4MNgB+GjnsCueo9zCLPEWJ+I9AGOB74NbiSt9gQ1HRsr2IG0UQaQD/yjtrQ0R0QSAh2EK9AbtIOILAS6NfLSXdT8Th2p+Uo6DpgnIhkmxIYFtfA73An8OLARtU5z8Rtj3qtd5i5qygSvBTK2YyCNPBdS+42vRCQReAu4yRhTEux4fCUiZwN5xpgVInJSsONpIxcwGrjBGPOtiDwJ3A7cE+ggwo4x5rSmXhORa4G3a5P5MhGxqJnfIT9Q8fmiqd9BRIYDfYHvRQRqShorRWS8MWZvAENsVnNtACAilwJnA6eG2gdrM3KA9HqP04DcIMXSZiLipiaxv2aMeTvY8bTSZGCGiEwH4oAkEfmXMeaXQY6rNXKAHGNM3TemN6lJ7gEViWWZd4FTAERkIBBDGE0+ZIxZa4xJNcb0Mcb0oWZHGR1Kib0lIjINuA2YYYw5FOx4WuE7YICI9BWRGGAmMD/IMbWK1PQI/g5sNMY8Hux4WssYc4cxJq12358JLAqzxE7tsbpLRAbVPnUqsCHQcYRlz70FLwEvicg6oAq4NIx6jpHiaSAW+KT228c3xphrghtSy4wxHhGZBXwEOIGXjDHrgxxWa00GLgHWisjq2ufuNMZ8EMSYotENwGu1nYRs4PJAB6BXqCqlVASKxLKMUkpFPU3uSikVgTS5K6VUBNLkrpRSEUiTu1JKRSBN7kopFYE0uSulVATS5K6UUhHo/wMGM0uLqohluAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "print('总体的均值为: ',data.mean())   #假装我们不知道\n",
    "print('空假设H0是：总体的均值是',data.mean() + 50,' \\n')\n",
    "\n",
    "sample_data=np.random.choice(a=data,size=80) #从中随机抽取80个做样本\n",
    "\n",
    "t_statistic,p_value=stats.ttest_1samp(a=sample_data,popmean=data.mean() + 50)\n",
    "print('从样本构造的t统计量 = ',t_statistic)    #t的绝对值越大，p值越小，可参见下面的t-分布图\n",
    "print('单样本双边检验的 p = ',p_value)\n",
    "\n",
    "#画一下样本容量为80时的t分布曲线\n",
    "df=80-1\n",
    "\n",
    "x = np.linspace(stats.t.ppf(0.00000001,df),stats.t.ppf(0.99999999,df),80)  #ppf函数是CDF的逆函数，用来求分位点\n",
    "plt.plot(x, t.pdf(x, df))  #pdf生成概率密度函数表\n",
    "plt.plot((t_statistic,t_statistic),(-0.01,0.4),'-.r')\n",
    "str_legend=('t distribution','calculated t')\n",
    "plt.legend(str_legend)\n",
    "\n",
    "\n",
    "x_025=stats.t.ppf(0.025,df) #注意：同样选取alpha=0.5，做双边test时要给左、右两侧各分配0.25\n",
    "plt.plot((x_025,x_025),(-0.01,0.4),'--g') #这根线是双边检验的左阈值\n",
    "x_975=stats.t.ppf(0.975,df)\n",
    "plt.plot((x_975,x_975),(-0.01,0.4),'--g') #这根线是双边检验的右阈值\n",
    "#print(x_975)\n",
    "\n",
    "rect1=plt.Rectangle((x_025,0),x_975*2,0.5,color='g',alpha=0.25)\n",
    "plt.gca().add_patch(rect1)\n",
    "\n",
    "rect2=plt.Rectangle((x_975,0),5,0.5,color='r',alpha=0.25)\n",
    "plt.gca().add_patch(rect2)\n",
    "rect3=plt.Rectangle((x_025,0),-5,0.5,color='r',alpha=0.25)\n",
    "plt.gca().add_patch(rect3)\n",
    "\n",
    "#plt.ylim(-0.02,0.5)\n",
    "plt.text(x_025,0.4,'accept region')   #图片上不支持中文\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p =  0.0941865267034535 \n",
      "\n",
      "\u001b[1;31m 接受 \u001b[0m H0\n"
     ]
    }
   ],
   "source": [
    "print('p = ',p_value,'\\n')\n",
    "alpha=0.05\n",
    "if p_value>=0.05:\n",
    "    print('\\033[1;31m 接受 \\033[0m H0')\n",
    "else:\n",
    "    print('\\033[1;31m 拒绝 \\033[0m H0,即总体的均值不等于',data.mean() + 50,'，此时错误拒绝H0的概率为',p_value,'小于显著性水平')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "不论是从计算得出的p值还是从自由度79的t-分布的图像中，都可以很容易的看出，t分布值处在接受域中，如果选择拒绝假设的话，有50%的可能性是错误的，因此我们没有充分的理由认为H0的假设不成立，选择接受假设。(这个重新运行过了，所以50%的概率值可能变了，不过不影响结论)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. 假设错误"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "总体的均值为:  8421.026851851851\n",
      "空假设H0是：总体的均值是 11421.026851851851  \n",
      "\n",
      "从样本构造的t统计量 =  -2.9088864530896843\n",
      "单样本双边检验的 p =  0.00470846938730825\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8VOX1+PHPuTPZSEjCFrawhVV22bcKbhXRIlardLFqXeqClp+lYt211lK1lvpV26q1VqtV6kqt4lI3FGRH9iUEAiFAQiAJCdlm7vP7YxIMkJBJcieznffrNSQzc+fOCc+9Z54597nPFWMMSimlIosV7ACUUko5T5O7UkpFIE3uSikVgTS5K6VUBNLkrpRSEUiTu1JKRSBN7kopFYE0uSulVATS5K6UUhHIHaw3bp+SYnp26hSst1cBVFJaiGmbGuwwopocKiQpUdsgEq3atu2gMaZDQ8sFLbn37NSJlX/9a7DeXgXQkuVvUnnZ94MdRlSLXfAmE8ZoG0QiOfPMbH+W07KMUkpFIE3uKmJ8+eGX7Nq+K2Drv+PqOygpLgnY+pVykiZ3FTG+/PBLsrc3/I3V6/E2af3z/j6PpOSkJr1WqZYWtJq7imx3X383efvyqKyo5JKrLuF7P/oeAMs/X85zjz6H1+slpW0Kj7/8OGWlZTxx/xNsXb8VQfjpL37K5PMns+KLFbww/wWqKqvo0r0Lcx+dS0JiAjMnzWTKBVNY/vly4uLiuPtPd3O44DBLPl7CN8u+4aUnX+KBPz9A1x5dj8Uzb848YuNi2b5xO4NHDeZnt/2MJ+5/gp1bd+L1eLnyF1cy6buTKC8rZ96ceezatotuGd04eOAgsx+cTf+h/Zk5aSZ/XfhXUtqmsOC5Bbz/7/cBuODyC7j0Z5eyP2c/c6+ay5BRQ9iwegMdOnbgoWcfIi4+Liht0BRVIuQkJlLucgU7lKgX7/WSXlpKTBOnZdfkrgLi9kduJzk1mYryCm646AbOOP8MjG147NeP8afX/kTnbp0pLiwG4MX/e5HE1ok8v+h5AI4UHaHoUBH/fPKfPPbPx0holcC//vIvFvxtAVfeeiUASa2TeH7R83zwxgc8+Zsn+d3ffseEcyYw/qzxTJ42uc6Y8vfn8+QbT+JyuXj20WcZMX4Ecx+ZS0lxCTdedCMjJ41k4T8X0jqlNS989AI7t+7k2guuPWk9W9dvZdHri3j6rafBwI0X38iwscNondKanF053PPEPcyZN4f7b76fL97/gnMvPjdA/8vOy0lMpHXHjvRMSUFEgh1O1DLGUFBURM6BA/QqaVop0K/kLiJTgT8BLuA5Y8y8epa7FPg3MNoYs7JJEamI8OYLb7L4g8UA5O/LZ++uvRQWFDJ0zFA6d+sMQHJqMgCrvlrFvU/ce+y1rVNas/R/S9mVuYtbLr0FAE+Vh4EjBh5b5qzpZwFw9vSzefqhp/2KafK0ybiqe6QrF69kycdLeO3Z1wCorKwkLzeP9SvXc8nVlwDQq38veg/ofdJ61q9cz6TvTiKhVQIAZ5x3BuuWr2PiuRPp3K0zfQb2AaD/kP7sz9nvV2yhotzl0sQeAkSEdikp5B882OR1NJjcRcQFPAWcC+QAK0RkoTFm0wnLtQZuBZY1ORoVEVbvyGHV6u089eZTxCfEM3vmbCorKjHGINSRNAwnJRNjDKMmjeKeJ+6p8z2OW97PPJSQkHDc+h94+gG69+5+0vs26BSLxMTGHPvdsiy83qbV94NJE3toaG47+HNAdQyQaYzJMsZUAq8CF9Wx3G+AR4DyZkWkwl5peSWtk1sTnxDP7h272bTG1w8YNGIQ3yz/hn179gEcK8uM+s4o3nrxrWOvP1J0hIGnD2TDqg3s3bUXgPKycvZk7Tm2zKfvfnrs56DTBwHQKrEVR0uP+hXj6DNG89Y/3jqWzLdv3A7AkFFD+Oy/nwGwa/susrZmnfTaoWOG8tWHX1FeVk7Z0TIWf7iYoWOG+vefo1QL8Se5dwX21LqfU/3YMSJyOtDNGPPuqVYkIteLyEoRWZlfVNToYJWDZs/23QJgbP8eeL1erjznSp75/TMMPN1XTkltl8ovH/4l995wL9ecfw0PznoQgCtmXcGRoiNcfd7VXHP+NaxZuobUdqnMfXQuv/nFb7hm6jXcfPHN7N6x+9h7HCk6wjVTr+GNF97g5ntuBuDM753Ja8+8xnUXXMfe7L2njPGnt/wUj8fDNedfw9XnXc3zj/vq/RddcRGFBYVcde5VPP+H5+nZtyeJrROPe22/wf0479LzuHHGjdx08U1ccPkF9B3U17H/v2hWWFjI088+69ey9z/8MI898QQA9z70EB9/+mm9y7797rts2rKl3uf/8re/8eIrrwAwZdo0Vq5e3eSYc/ft49IrrvD79YEiDX0NFZEfAOcZY66tvn8FMMYYc0v1fQv4BLjKGLNLRD4D5jRUcx/Vv7/RM1SDqCaxz5/v6Go/PryKjVsWM+yHgfngAI4bteKU/GKbFxdXUlLmwfZ6ccXEYpXsY91f7+TlT186rtwS6lblrsL9+WJ+Mb7xbbA5JYXT+vQJQFT+2ZWdzYWXXcaGZQ1Xd+9/+GGSkpKYc+utDS571Q03cOHUqVw6Y8ZJz3k8HtzubyvUU6ZN47GHHmLUiBGOx9xYmzMzOe2EjrCceeYqY8yohl7rT889B+hW6346kFvrfmtgMPCZiOwCxgELRaTBN1eR56Hsl/iHJ7yOpVdUGZ74oILMA17ipZJvnvoVKx+/mcVP/YYuF9wErvAaVPbSuvBrgxp33HcfO3buZPjEifzq7rtPev63jz5Kv9NPZ9J3v8vW7duPPX7VDTfw+ttvH1vHwNGjGTp+PHPuuosly5ax8L33+NU99zB84kR2ZGUxZdo0Zs+dy6jJk/nTn/983LcAgJdefZXhEycyeOxYlq/0/V+euMzgsWPZlZ19Usy7srMZPHYsAOXl5Vx9440MGTeO0ydN4tMvvgDghZdf5vs//jFTL76YvsOHc/s9dR9bag5/ttoVQF8R6QXsBWYCP6p50hhTBLSvue9vz12ppnr1y1cdW5cxhr99VkFOgc3/mxbH0O6JcMkzACzf4eHpjyr455eVXHVGbNQdaHzgwx1sOuDsGbkDOyZx33dPHoFUY94DD7Bh82bWfvXVSc+tWrOGV994g7VffYXH42HEd77DyNNPP26ZQ4cO8dZ//sOWVasQEQoLC0lNTWX6tGkn9dwrKytZ+fnngC9x13a0rIy1X33FF199xc9uvvmUvfITY96V/e2JdE9Vl2vWf/01W7Zt47szZrCtuuSzdv161ixeTFxcHP1HjuSWn/+cbunp9b5PYzXYczfGeIBZwAfAZmCBMWajiDwoItMdi0SpIHh/bRXLd3i5dGwMQ7sf39cZ09vNhafH8PlmD59u8gQpQlVj8ZIlXHzhhbRq1Yrk5GSmT5t20jLJycnEx8dz7axZvLlwIa1atap3fZdfckm9z/3w0ksBOGPiRIqPHKGwsLBJMX+5dClXzJwJwIB+/ejRrRvbMjMBOHvyZFJSUoiPj2dg//5k79lzqlU1ml/fN40x7wHvnfDYvfUsO6X5YSkVeOt2e/j3sirG9HYxbXjdNfXvj45hT4HNy19V0rWtRf/O0XPm5ql62MHS0Lcnt9vN8k8/5X+ffcbr77zDk888wyfv1j3OI/EUif/E9xER3G43tm0fe6y8vOGBgac6phkXG3vsd5fLhcfjbAdC55ZRUam0wvCXjytIb2dxzZS4epOGZQk/PzuODq2Fpz6soKKqaaeCK/+0TkriSD1nZJ4xcSJvvfsuZWVlHDlyhP+8//5Jy5SUlFBUXMy0887jj7/7Hd+sX9/geuvy2ptvAr6ed0pyMikpKfTs3p3Va9cCsHrtWnZWl18aivnlBQsA2LZ9O7tzcujft2VGVoXXkSLlnNtuC8hq/9rvNlav+ygg63bSl1s8HK2E26fEEhdz6t5gqzjhqslxzFtYzteZHiafFtojZ24bfxsx74d+G9SlXbt2TBw7lsFjx3L+uefy6EMPHXtuxPDhXP797zNswgTSOnRgdB2jWY6UlHDRzJmUV1RgjOHx6lr6zEsv5bpbbuGJv/yF1198scE44uPiOH3SJKqqqnj+qacAuOSii3jxX/9i0JgxjB01in7Vo4pOjPnm6647tp6brr2WG2bPZsi4cbjdbl7485+Ji2uZuYYaHAoZKDoUMnKF+sU6bGP49b/KaJ0g3H1xQsMvwPf1+p5/l+GyhPsviQ/5g6tNvVhHsIdCquMFeiikikRLlvhuDvvPwSV86d3p+HqdtDHHy4Fiw9mD/e+BiwhnD44h+6DNjgN2wy8IoiV7Qr8NVOBpco9WCxb4bg77Q84CXvWscXy9Tvp4g4fkBBiV0biDo+P7ukmIhY83VAUoMmcs2Bj6baACT2vu0eqBB4IdQVDkF9usy/Zy4YgYYlyNK63ExwiT+rn5ZJOHHx41pLQK7dKMim7ac49WKSm+W5T5dJMHEThzYNP6NWcNjsFrw+ebQ7v3rpQm92i1aJHvFkUqPYbPN1cxopeLtklN2/Q7p1oMSrf4bJMHr63DIlXo0uQeraIwuS/L9FBaAWcPat5QxrMHxXCo1LBmV/jN1a6ihyZ35aiXBtzJPTGheVm5TzZ66NJGGNCleZv98B4u2iUJn2wMzdLMnd8J3TZw2gsvv8ysX/6ySa+tPdnYqdafu29fo9Zbe+KwEx9/JQCDGOqjyV05qlt8Gh2t1sEO4yQFR2x25ttM6u9u/hVuLGFCPzebc21KykOvNJOWGJptEI6aktzrs2v3bl75978dWZc/NLkrR72W9wkfe7YFO4yTfLPbV0IZ3sOZAWKn93BhDKzbHXqlmU92hmYb+OvFV15h6PjxDJswgSuqz/b8z/vvM/bMMzl90iTOmT6dA3l5J73uQF4eF//oRwybMIFhEyawZNmyk3rRjz3xxEkzQAI8OG8eoydPZvDYsVx/660YY3j97bdZuWYNP772WoZPnEhZWRmr1qxh8vnnM/KMMzhvxgz27fddI3fVmjXH3vepei42csd997F46VKGT5zIH5980on/qlPSoZDKUX/OXUix9yBnBDuQE6zd5aVjstA51Znhiz3TLJIThLXZHib0C63daOHWhVjeg9Q5s19j1THz4kmmToWaC2ZMmwY//rHvVlAAJ16R6L33Tn59LRs3b+a3jz3GVx99RPt27Th06BAAk8aN4+tPPkFEeO4f/+CR+fP5wwlJ+tbbb2fyxIm89coreL1eSkpKOOznbI6zrr+ee++4A4ArrruOdxct4tIZM3jymWeOXbijqqqKW371K9559VU6tG/Pa2+8wV0PPsjzTz/N1TfdxP89+iiTJ02qcx568E0N/NgTT/BuC/XeQ2urVCoAKqoMm3K9nDWw+SWZGpYIw3q4WJnlweM1uBs5Zl7V7ZPPP+fSGTNo364dAG3btgUgJzeXy6+6in0HDlBZWUmvHj3qfO2L1VOauFwuUlJS/E7uny5ezCPz53O0rIxDhw8z6LTT+N755x+3zNbt29mweTPnXuS7hLTX66Vzx44UFRVRWFTE5EmTALhi5kze/yj4c/toclcRb2OOF48Xhvd0dnMf3sPF4i0etu+3Oa1rhE4F3EBP+5TLt2vX6NcbY6jrY/KWX/2K22bNYvq0aXy2eDH3/+53fq3Pn2l6y8vLuem221j5+ed0S0/n/ocfrnM5YwyDBgxg6f/+d9zjhYWFITnXkNbcVcRbm+0lIRb6dXJ2cx+U7sJtwdpsvZCHU86eMoUFb71FQUEBwLGyTFFxMV07dwbgH9UXsj7ptZMn8+fnngN8veri4mI6pqWRl59PQUEBFRUVvFvH8N+aRN6+XTtKSkp4/Z13jj1Xezrf/n37kn/wIEurr8pUVVXFxs2bSU1NJSU5mS+XLgU4NsXviRo77XBzaXJXEc02hnW7vQxOdzleOomPEQZ0dfFNdugdVA1Xg047jbvmzGHytGkMmzCB2+68E4D7f/1rfnDllYw844xjJZsT/emRR/h08WKGjBvHyDPOYOPmzcTExHDv3LmMOfNMzpk+nQH9+p30utTUVK678koGjx3LeRdffNxUwlf9+MfcMHs2wydOxOv18vpLLzH3vvsYNmECwydOZEl1ov/7009z8y9/yfCJE+u9QMfQwYNxu90MmzChRQ6o6pS/0apmGlGHpyA4WFXE8tXv0uryHzu63qbameflgTfLue6sWCb2c34e9o83VPHPLyuZNzOBTqmh0VcqKi8i5p13mTa28W2gU/6GFp3yVzVegOaWaR+TQqr4N0d6S1ib7UWAod0Cc3hpWHfXsfcJFSnxodUGKjg0uUerAE0/8ML+RfzXs9nx9TbV2mwvvTtatE4IzAGvDskW6W0lpOruizJDqw1UcGhyj1YBTO7ve0MjsRwusck+aDO8Z2BHsgzv4WbbPpvSitA4W3VRZvPaIFilWnW85raDJvdoNX++7xbBnD4rtT7DeriwDazfEzqlmaaK93opKCrSBB9kxhgKioqI9zZ9m9Jx7ipifbPbS7skoWubwI5B7p1m0Toevsn2MK5PeO9S6aWl5Bw4QP7Bg8EOJerFe72kl5Y2+fXhvSWqpnvtNd/Pyy8PbhwBYtuGLblexvR27qzU+liWcFpXF5v32r6TcELwhBZ/xRhDrxYci60CR8sy0WrpUt8tQu05ZFNWCf07t8yZowO6uCg8asgv1nKGCg3ac1eOem/IPL5e+U7DCwbYllzfKef9mzl3u79qPkS27POSlhLcPtO8c+YR+0bw20AFl/bclaNaueKJF+dPFmqsrbleOrQW2jXxcnqN1aWN0DoetubaDS8cYPHu0GgDFVya3JWjnt77Nm961gU1BtsYtu330r9Ly03mJSL06+xi677gj5h5e0vw20AFnyZ35agF+Z/xiTczqDHkHjaUlEP/zi27eQ/o4uLgEUPBkeD23j/bFfw2UMGnyV1FnC25vt7zgBbsucO3HyZbQqD3rpQmdxVxtu3z0jZRaN+6ZYckpre1aBUbGnV3pTS5q4hijGFLrk3/LlaLjze3rNCpuyulyV1FlP1FhuIy02Lj20/Uv4uLA0WGwlLtvavg0nHu0SpA88p8Nnw+S5a/SWVA1t6wrdX19pYcKVPbt3V3m3F9gtN3mj91PrEL3gzKe6vQoT13FVG27vOSnCB0SgnOFAA92lvEx8C2XC3NqODS5B6tXnvt2/llHPTYntd4pWq14+v1hzGGrbk2A4JQb6/hsoS+nVxBHTHz2obgtYEKHVqWiVYbNwZkte8WLKXYPsilAVn7qeUfMRwqDV69vUb/LhavL/NSXGZIDtBFQk5lac5SLFtndYx2mtyj1YMPBjsCxwW73l7D9+FSxbZ9XkZl6C6mgkPLMipibNtnkxTvm+clmHp1sIh1o0MiVVD5ldxFZKqIbBWRTBG5o47nbxCR9SKyVkS+FJGBzoeqHPXss75bBNmR56V3RxdWkOdTd7uEXh0sdhzQ4ZAqeBpM7iLiAp4CzgcGAj+sI3m/YowZYowZDjwCPO54pMpZGzcGpO6eYMURR8uXRUorDLmHDb3TQuPLaO+OLrIP2lR6Wn5+9zh3cNpAhRZ/9oQxQKYxJssYUwm8ClxUewFjTHGtu4mAXrEgSr0/9Pf8Ie6ihhd0WFaerwTSp1NoJLXeHS28Nuw+2PK999+fE5w2UKHFn+TeFdhT635O9WPHEZGbRWQHvp77rXWtSESuF5GVIrIyv6ioKfEqVacdB2wEX707FNR8g9DSjAoWf/aEugqYJ/XMjTFPGWN6A3OBu+takTHmGWPMKGPMqA4pKY2LVIWF32S/yN+rlrf4+2YesOnaVkiIDY3rl6YmWrRLEjIPtPxB1Re/CU4bqNDiT3LPAbrVup8O5J5i+VeBGc0JSoWv/x1ezSo7p0Xf0zaGrAO+g6mhpE8nix15Ld9zX72v5dtAhR5/kvsKoK+I9BKRWGAmsLD2AiLSt9bdC4DtzoWo1KntLzQcrYQ+HUOjJFOjd5qLQyWGwyVamlEtr8EzLIwxHhGZBXwAuIDnjTEbReRBYKUxZiEwS0TOAaqAw8CVgQxaqdp2VJc+Qq3n3rv6wyYzz2Z0C13LVakafp0+Z4x5D3jvhMfurfX7LxyOSym/7ThgkxALnVJDo95eo0d7C7fLF9/ojGBHo6KNnhsdrZKTA7LadjHJQElA1l2fHXk2vdOCf/LSidwuoWd7i6wWPqiaHJeM1cJtoEKPJvdoFaC5Zd4Y9GCLzudeVmnIOWQzYkRMC71j4/TuaPHJRg8er8HtapkPnwfPfFDnc1c6t4wKbzvzbYz5tr4danp3dFHlhT0FelBVtazQ3CNU4AVobplfZz3Ln6uWOL7e+tQcTM1IC62DqTVqPnRackjks6tatg1UaNLkHq2Kinw3hy0t3shGe7/j663PjgM2nVKFpPjQqrfXaJsopLaSYx9CLWFjfsu2gQpNWnOPVnPmBDuCZjPGsOOAl2E9QnczFhF6d7TI1GkIVAvTnrsKW/nFhiPloVtvr9Gno4v8YkNxmc6np1pOaO8VKnAee8x3C2M1veFQmea3Psfq7kGYZ0ZFr9DeK1Tg5OT4bg5Lj+tAB0lyfL11ycrzEuuGrm1DezPu2cHCEshqoYOqHRJbrg1U6ArdYqUKS/887a4WG+eelWfTq4OFywrNg6k1Yt1Ct3bWsTnnA+2u79xF7D4d5x7tQrvLo1Q9PF7D7oM2vUJ0COSJeqVZ7MyzsY3W3VXL0OSuHDU780nmV34R8PfZU2DjsSEjxOvtNTLSLI5WQl5R4JP7k8tbpg1UaNOyjHLU2pJMis3BgL9PTf06fJK77xtGVp5Np9TAxpx5KBOrBdpAhbbw2DOUOkFWnk1yArRLCu16e40uqUKcmxaruyulyV2FpZ35XnqluZAQmwmyPpYl9OxgtdiIGaU0uauwc7TCsO+wCZuSTI2MNBe7D9p4vHpQVQWe1tyjVXp6QFbbr1U6B0o8AVl3jV0HbQzhU2+vkZFm4bF9B4MDOconPTkd18HAtoEKfZrco1WA5pZ5pt8clhQGdpx7Td26V4fwGAZZo+bDKCsvsMl9zoQ5xOboOPdoF15dH6WAnXk2HZNDdybI+rRNEpITROvuqkVoco9WAZpb5vptj/H7yk8cX29tvp5v+G26IkJGmsXO/MCOmHlsSeDbQIU+LctEq5SUgKx229Ecik1hQNYNcLjE5nCpCdmLczQkI83im2wvRysMreIC880jpzgHK4BtoMKDJvdodd11wY6gSbLyw+vkpRNlpFkYfAeFB3YNzw8oFR7Ccw9RUWtnno3Lgu7tw3PT7dmh5kxVPZlJBVZ47iGq+e6913cLM1l5XtLbWsS6w+tgao2keKFjirBTD6qqANOyTLQqLg7Iaocn9SG3JCCrxjaGnfk24/qE92abkWaxJTdwyb1P2z5YOrVM1NOeu3LU/D6zmB17RkDWvb/QUFYZvvX2GhlpLg6XGg6XBCbBzxoTuDZQ4SO89xIVVXbWnLwUpiNlatQM46w5OKxUIGhyV476yebf8kDlhwFZ9448m/gY3wyL4ax7OwuXFbjL7v12ceDaQIWP8C5eqpCTU5FPsQlM0b3msnpWiF9WryGxbqF7O4usAF0wO780HytAbaDCh/bcVVio9Bj2FNhkdAzvkkwN35mqNratM0SqwNDkrsLC7oM23jC6rF5DMtIsyqtgX6EmdxUYkbGnqIi3I8wuq9eQmm8gO/RkJhUgkbGnqMYbNMh3c9j45EEMsjo5vt6deV7aJgptEiNjk+2YIrSKDcxB1UEdAtMGKrzoAdVoFaC5ZX6XcR1LDjo/n/uOPJuMjpGR2AEsEXqlWWQdcD65XzfyOmJ36Hzu0S5y9hYVsYrLDPnFhowOkbW5ZqS5yDlkU1GldXflvMjaW5T/AjS3zCUb7+XOivccXWfNyUuRMlKmRkaahW0g+6Czvfd7P3W+DVT40bJMtApAvR2goKqYYsodXWdWno0I9IywnnvvYwdVbfp1du6Dq7iiGMvhNlDhR5N7tLr88mBH4LcdeTZd2wjxMeF98tKJkhOE9q2l+mSmmGCHoyJMZHWFVMQxxrAzzxu2V15qSEaapddUVQHhV3IXkakislVEMkXkjjqev01ENonIOhH5n4j0cD5U5ajZs323EHegyFBaAb0jaKRMbRlpLgpKDEVH9aCqclaDe4yIuICngPOBgcAPRWTgCYutAUYZY4YCrwOPOB2oCg9ntxnBSCvdsfVlHTt5KTJ77r1rZoh08GSmEZ2dbQMVnvzpDo0BMo0xWcaYSuBV4KLaCxhjPjXGHK2++zWgW1aUuqfHT7k6Zoxj68vK8xLnhq5tIqveXqN7ewtLnD2Z6afDnG0DFZ78Se5dgT217udUP1afa4D3mxOUUjWy8mx6RsBMkPWJixG6BXCGSBW9/Enude1VdRYIReQnwCjg0Xqev15EVorIyvyiIv+jVGHj/HVz+WXFO46sq8pr2H3QjtiSTI2MNIusfBvbOFN3n/uxc22gwpc/yT0H6FbrfjqQe+JCInIOcBcw3RhTUdeKjDHPGGNGGWNGdUhJaUq8KsSV2RVU4EwvdPdBG49NRE07UJeMNIuySth32JnkXuFxrg1U+PJnr1kB9BWRXiISC8wEFtZeQEROB/6KL7HnOR+mikaZ+3116D4Rntz7dPJ9M8nU0oxyUIN7jTHGA8wCPgA2AwuMMRtF5EERmV692KNAEvBvEVkrIgvrWZ1Sftt+wEv71pEzE2R9OqUIiXGQGYBJxFT08usMVWPMe8B7Jzx2b63fz3E4LhXljDFk7rcZ0CWyEzuAiNCnk4vM/dpzV87R6Qei1fjxAVnthe3Gs6t0fbPXU1BiKDxqjpUsIl3fjhbfZHspKTckxTdvZND49PG4CprfBiq8aXKPVgGaW2ZOt8tZsi+m2fO515Qo+naK/J471NTdq9hxwMuwHs3bLS8ffDmxm3SummgXHXuOCjuZ+30nL6W3jY5NtFcH38lMWndXTomOPUedLEBzy0xZO5tZFc2/ClDmAZuMNAtXhJ68dKK4GKF7e8uRuvvsRc60gQpvWpaJVlOnBjuCelVU+U5k2BjGAAAV1UlEQVReuuD06Cot9Olo8cUWD17bRM2Hmgoc7blHq6lTQzbBZ+XZ2Cbyx7efqE8nF5Ue2FOgpRnVfNG196hvFRX5biGo5mSe3hF2Wb2G9K3+MNO6u3KCJvdodd99vlsIytxv06WNNHtIYLhpmyS0SRQd764coTV35ajLOkwh6+jaJr/eNobMA15G9oq+TVNE6NPRYvv+5vXcp/ScgvtQ09tARYbo24NUQN3UdQZL9tpNHue+v9B35aU+UTK+/UR9OrlYkeXlcIlNm6Sm/R/MGDCD2HVa2ol20bkHqYA56i2n3FQ1+fU19fY+UVZvr9HHgbp7uad5baAigyZ35ahp6+9gTuV/mvz6zP02iXHQKTW66u01erS3cLuaN0PkHR83rw1UZNDkrkJK5gEvvTu6sCQ6k7vbJWR0sHTEjGo2Te4qZJSUG3IPm6gb336iPp1c7Mq3qahy5uIdKjpF916kQsrWfb5SxIAu0Vlvr9G/s4XXhh3ae1fNoMldhYwtuV5i3dArLbo3y36dXYjAln063l01nQ6FjFYBmnrgqk5T2X50VZNeuyXXpk9HixhXdNbbayTECj3bW2zJbVpyn9pnKq5DTWsDFTk0uUerACb3JbuPNnqce0m5IafA5uLR0TVZWH0GdHHx0foqKqoMcTGN+7Cb2mcqsauPBigyFS6i+/tvNAvQ3DIHq4ooNGWNft3WfV4MWm+vMaCLhceGHXmNr7sXlTetDVRk0eQerQI0t8ylG+/j7sr3G/26LbleYlxab69xrO7ehNLMfZ81rQ1UZNGyTLS67LJgR3Ccrbk2fTppvb1GTd19axPr7kppNylaTZjgu4WAknLDngJbSzIn6N/FxY4DNpUeHe+uGk+Te7Tavdt3CwHbtN5ep2N1dx3vrppAk3u0evxx3y0E1NTbM7Tefpx+nZped1dKa+7KUTd2mc7WzOWNes3WfTq+vS6t4oQeTRjvPr3/dNyHGtcGKvJoV0k56vK0szjH3c/v5UsrfBfD1pJM3QZ0sdiR17i6+1m9GtcGKjJpcleO2lOexwH7iN/La7391AZ0ceHxNq7unlfauDZQkUmTu3LUFVse5jdVH/m9vI5vP7W+Tai7P7y4cW2gIpPuUSqoauaTiXVrvb0uiU2suyulyV0FTdFRQ/ZBm9O6aknmVAZ2dZF5wKasUse7K/9pcldBs2GPB4Ch3TW5n8qQbi68Nmzaq7135T9N7ipo1u3xkpwgdG+vm+Gp9O1kER8D63drclf+03Hu0SpAc8v8Mv0yNm9f2uBytm3YsMfL8B7uqL1eqr/cLmFQuot1e7wYY5AG/r8uG3QZ7i8bbgMV2TS5R6sAzSvzvfYTaJe1v8H53LPybEorYIiWZPwypJuLVTu95B42dG176uQ+odsEYl37WygyFar0+3C0CtDcMluP7ibbPtzgcuv3eBGBwema3P1R8yG4bk/DpZndRf61gYpsmtyjVYDmlvn5tsd5tOrTBpdbt9tL7zSLpHgtyfijXZJF1zbC+t2eBpd9fKl/baAim5ZlotW11wbtrYvLDLvybWboJfUaZUh3Nx+vr6K8yhDfyEvvqeijPfdoNXiw7xYEG/b4phwY2k1LMo0xtLsLjw2bdUik8oMm92i1YYPvFgTrdntIToAeHXTza4x+1UMi1+mQSOUHv/YuEZkqIltFJFNE7qjj+TNEZLWIeETkUufDVI577jnfrYXZtmFDjpfB3XQIZGO5XcJpXV2s2+0bEqnUqTSY3EXEBTwFnA8MBH4oIgNPWGw3cBXwitMBqvByd48ruNI9qt7nd+bblJRrSaaphnZ3UVBi2FdYf3K/Yuip20BFB38OqI4BMo0xWQAi8ipwEbCpZgFjzK7q5/R6YFHunDYjaeXKrnec+7rdXgQYpEMgm2RI9Yfiut1eurSpu282sstIYl3ZLRmWCkH+lGW6Antq3c+pfkypk6wtyWSbnV//89leMjpatE7QkkxTtG9t0aWNsDa7/iGRmYdO3QYqOviT3OvaC5tU8BOR60VkpYiszC8qasoqVIibnfkkT1QtrvO5vCKb7IM2ozJ0BG5zjM5ws3WfTeHRur8oP7m8/jZQ0cOf5J4DdKt1Px3IbcqbGWOeMcaMMsaM6pCS0pRVqDC2fIevtzk6Q0syzTGmtxtjYGWWjppR9fMnua8A+opILxGJBWYCCwMblopEy3d46d3Ron1rHQLZHF3b+kozK3Y0fLaqil4N7mXGGA8wC/gA2AwsMMZsFJEHRWQ6gIiMFpEc4AfAX0VkYyCDVuFnf6HN7gKbMVqSccSY3m627bM5XKpjGFTd/NrTjDHvAe+d8Ni9tX5fga9co1SdjpVkemtJxgljert5e2UVK7O8nDtEvwmpk2k3KloFaG6Zh3tdy/pNn5/0+IosL307WbRN0kTkhC5tLNLb+koz5w45fo6ea0dcS8z/Tm4DFV10T4tWAZpbZkLKYIa4Oh/3WO5hmz0FNmN6a1/CSWN6u9m23+ZQyfGlmcFpJ7eBij6a3KNVgOaWWVK0gfXefcc9tmKHBwFG9dKSjJNGV39Yrjhh1MyGvJPbQEUfTe7RKkBzy9y58zn+6jn+Em/Ld3jo29mijZZkHNU51aJbO+ukUTPPrT65DVT00e/J0eq221rkbfYestl72PCTSTp3eyCM7e3i9eVVFByxaadDTFUtujVEq+7dfbcAW1ZTktETlwKipjSzTMe8qxNoco9WS5b4bgHktQ2Lt3gY1M1Faivd1AKhY4pFn44WX2z26DTA6ji6x0WrBQt8twBam+3lcKnhrIFa/QukMwe52V9k2LxXT2hS39Lkrhw1v88sbo35DgCfbKyibaIwrIeWZAJpdIabpHj4ZFMVALPGfNsGKnppcleOGp7Uh35WB/YX2mzMsZk80I3L0ul9AynWLUzqH8PqnV4Ol9r0aetrAxXdNLkrR318eBUrvLv5bHMVLgsmD9CSTEs4c6Ab28AXmz2syvW1gYpuuucpRz2U/RKFVQep2uLh9J4uUhO1/9ASOqZYDE538dlmD/bBl3B5DvKLYAelgkr3POW4YhNLaQWcPUjHtrekswa5OVxqKCnXUTNKk7sKgMN2HJ1ThQFddPNqScN6uGibKBSWanJXmtyVw0q9hjLj5syBMYjogdSW5LKEyQPdlFYYKozu2tFOtwDlqNwyg4VhYn89nBMMkwe4EYECOz7Yoagg0z0wWgVgbpnNR2wovZkfxO0mMU577cGQmmhxYY9fsHR7FbtKbXrqAe2opS0frQIwt8zj2ytp50rnJ/Fljq5XNc6Px2XQynTlTzsqgx2KCiJN7tHK4bll1hZ6+SjPy/hOq1hjtju2XtV4mwq+pl/8f3k718v2Ep2SIFppco9WDs8t84ftVbSNgU1Vb/KqZ41j61WNt2DjAnZb75Hogj9u1957tNKae7R64AHHVrXskJfFBV7u6h/Ly0WOrVY1g1sMP+kZwxM7qthQ7GVwss7vE2205x6tUlJ8t2YyxvCH7ZWkxQk/6a59hVByTc8YUmLgj9urgh2KCgJN7tFq0SLfrZk+zvey/LDNrIwYElw6QiaUpMQI1/eM4X/5XpYUeBt+gYoomtyjlQPJvajKcNfGSgYkCTO7aa89FP2sZww9WwlzN1Rw1KNnrkYTTe6qyR7YXElBpeGxIXHEVk/r+9KAO7kn5twgRxbd7vzOt22Q4BIeGRxHTpnh99v04Go00eSumuTjPA9v5nq4OSOGwSnfHqzrFp9GR6t1ECNTaYnHt8GYti6u6uHmH7s9LNXyTNTQ5K4arajKcGd1OWZW7+Nnfnwt7xM+9mwLUmQK4JOdJ7fB7f1i6dFKuF3LM1FDk7tqtLrKMTX+nLuQt70bghSZAli49eQ2SHAJj2p5JqpocleN8sqeKt7M9XDTCeUYFfpql2feyfUEOxwVYJrcld/+l+fh7o2VTGnv4tbeeiGOcDS3Xyxj2ljMWV+hwyMjnCZ35Ze1hV5mfVPBoGSLp4bHEaMXvQ5L8S7h2RHx9EoUfr6mnC1HdO6ZSKXJXTVoV6nNNavLaR8rPD8ynkS3JvZwlhIj/H1kPK1cwlUry9lXrgk+EumZJ9HKz7llskptrlpZjm3gH6Pi6dDAPO2vD3qA5avfdSJC1UQPTHmAmHdO3QZdEyz+PjKOy5aVc8WKcv4+Mp5urbSvF0m0NaOVH3PLLD7oZcbSMko8hr+PjCfDjws/tI9JIVUSnIpSNUFKvH9tMDDZxXMj4zlQYZjxdRkrDmsNPpJoco9Wp5h+wBjDP7KruGpVOZ3jhXfGJzA81b+RMS/sX8R/PZudjFQ10qJM/9tgXFsXb49LINkt/Gh5OQtydJKxSKHJPVrVk9wLKw1zN1Ry3+ZKzmzv4o1xCY36uv7C/kW879XkHkyLMhvXBr2TLN4en8DYtha3b6jknk0VHNETncKe1tyj1fz5x921jWFBjodHtlVSWAU3Z8RwW98YXKIHT6NBzUHW322t5O/ZHhYd8HJn/1hmdHYhug2EJe25RzljDF8f8nLx1+XcsbGSPkkW/52YwK/6xWpijzIxlnDvaXG8PT6eLvHC/1tXwWXLy1l12Isx2pMPN9pzj1IV/3qVDUU2dw2ewZYSQ4c4Yf7QOC7SnlrUG5bi4q1x8fx7r4ffb63kkmXlDE62+Gl3N9/r7NZ5+8OEXz13EZkqIltFJFNE7qjj+TgRea36+WUi0tPpQFXzFVYaFu7zcNu6CtYt+oqqJUtxWcLvB8fyxRkJzOji1sSuALBEuDw9hi8mt+I3A2OptA23b6hk3GdH+dX6Cv6730NRlfbmQ1mDPXcRcQFPAecCOcAKEVlojNlUa7FrgMPGmD4iMhP4PXB5IAJW/vHYhqxSw8YjNhuLvawutFlbaGMDqTFwU6zQMU54d3y8own9vSHz+HrlO46tTzXevHPmEfuGM22Q5Bau6B7DT7q5WX7Y5uU9VSw64OHfez24BEakWpye4mJQssWgZIteiaLlvBDhT1lmDJBpjMkCEJFXgYuA2sn9IuD+6t9fB54UETFaqHOUMYajXij1GEq8cLjScKjKcKjSUFBp2FtmyCkz5JTZ7CkzVFafeBhnwcDWFrN6xzClg4thKRau/1R/aXN4R2zliideYtB5B4Mn3h1PrDg794+IMLati7FtXXhsw5oim8/zfRdGfyG7isrqPT3egvQEIT3BolsroUu80C5WaFt9axMjJLohySUkuNBvigHkT3LvCuypdT8HGFvfMsYYj4gUAe2Ag04EWdvHeR7ebsEZ7fz5dKr5CKtZtvZPY6p/Vj9gABuwTfUNg7f6d6+BKhuqDHiMocqGChsqvIYKG456Tx1PmxhIT7Dol2RxdgeLAa2FQckueicK7haaC+bpvW+T5VnHhXy/Rd5PneztLW/j9qxjQoDawG0Jo9u4GN3GxRygyjZklho2FHnZcsQmp9zXyVhT5KXoFMPmLSDB5et8xLmEOAtiLHCLEGuBW8AlYNX8RLCq71uA4PtHqn+X6t+hjp9+bP4t+TFzeXoM32kf2FlV/Unudf3NJ+YYf5ZBRK4Hrgfo3r497N/vx9sfr6Aghk2HW3ZGwlM1+okbT70bV/UyFtUbZvWG6BJwYbAEYoDWAm4xxFgQ44Z4yxBnQaxAK8uQaEGiy5BkGVLdhnZuQ1u3oa3L0KqubaW0+naiyuq+dRPa4FQW7P2QI1Y5o0qcXa/y34c7PiTBKne8besTA5wGnOYG2hz/XKkXDnmEQ17hkEco9AgltlDqFUpsKLPF14ExQrkNHuPr3FTZgseAF1/HpwrBa77tMNn4bjWdJzjFT9NwJ62lSwyH40rAE9hOqjRUORGR8cD9xpjzqu//GsAY87tay3xQvcxSEXED+4EOpyrLjBo1yqxcudKBP0E1yZQpvp+ffebsal/wrfezq5xdr/KftkFkE5FVxphRDS3nz2iZFUBfEeklIrHATGDhCcssBK6s/v1S4BOttyulVPA0WJaprqHPAj4AXMDzxpiNIvIgsNIYsxD4G/CSiGQCh/B9ACillAoSv05iMsa8B7x3wmP31vq9HPiBs6EppZRqqgZr7oGiNXellGo8J2vuSimlwowmd6WUikCa3JVSKgJpcldKqQikyV0ppSKQJnellIpAmtyVUioCaXJXSqkIpMldKaUiUNDOUBWRfCDb4dW2JwBzyLcw/RuCL9zjB/0bQkUg/oYexpgODS0UtOQeCCKy0p/TckOZ/g3BF+7xg/4NoSKYf4OWZZRSKgJpcldKqQgUacn9mWAH4AD9G4Iv3OMH/RtCRdD+hoiquSullPKJtJ67UkopIjS5i8gtIrJVRDaKyCPBjqepRGSOiBgRaR/sWBpDRB4VkS0isk5E3hKR1GDH5C8RmVq97WSKyB3BjqexRKSbiHwqIpurt/9fBDumphARl4isEZF3gx1LU4hIqoi8Xr0fbBaR8S0dQ8QldxE5E7gIGGqMGQQ8FuSQmkREugHnAruDHUsTfAQMNsYMBbYBvw5yPH4RERfwFHA+MBD4oYgMDG5UjeYBfmmMOQ0YB9wchn8DwC+AzcEOohn+BCwyxgwAhhGEvyXikjtwIzDPGFMBYIzJC3I8TfVH4HYg7A6KGGM+NMZ4qu9+DaQHM55GGANkGmOyjDGVwKv4Ogphwxizzxizuvr3I/iSStfgRtU4IpIOXAA8F+xYmkJEkoEzgL8BGGMqjTGFLR1HJCb3fsB3RGSZiHwuIqODHVBjich0YK8x5ptgx+KAnwHvBzsIP3UF9tS6n0OYJcbaRKQncDqwLLiRNNp8fB0bO9iBNFEGkA/8vbq09JyIJLZ0EO6WfkMniMjHQKc6nroL39/UBt9X0tHAAhHJMCE2LKiBv+FO4LstG1HjnCp+Y8w71cvcha9M8HJLxtYMUsdjIbXd+EtEkoA3gNnGmOJgx+MvEbkQyDPGrBKRKcGOp4ncwAjgFmPMMhH5E3AHcE9LBxF2jDHn1PeciNwIvFmdzJeLiI1vfof8lorPH/X9DSIyBOgFfCMi4CtprBaRMcaY/S0Y4imdqg0ARORK4ELg7FD7YD2FHKBbrfvpQG6QYmkyEYnBl9hfNsa8Gex4GmkiMF1EpgHxQLKI/NMY85Mgx9UYOUCOMabmG9Pr+JJ7i4rEsszbwFkAItIPiCWMJh8yxqw3xqQZY3oaY3ri21BGhFJib4iITAXmAtONMUeDHU8jrAD6ikgvEYkFZgILgxxTo4ivR/A3YLMx5vFgx9NYxphfG2PSq7f9mcAnYZbYqd5X94hI/+qHzgY2tXQcYdlzb8DzwPMisgGoBK4Mo55jpHgSiAM+qv728bUx5obghtQwY4xHRGYBHwAu4HljzMYgh9VYE4ErgPUisrb6sTuNMe8FMaZodAvwcnUnIQu4uqUD0DNUlVIqAkViWUYppaKeJnellIpAmtyVUioCaXJXSqkIpMldKaUikCZ3pZSKQJrclVIqAmlyV0qpCPT/AT+QG4hd08CRAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "print('总体的均值为: ',data.mean())   #假装我们不知道\n",
    "print('空假设H0是：总体的均值是',data.mean() + 5000,' \\n')\n",
    "\n",
    "sample_data=np.random.choice(a=data,size=80) #从中随机抽取80个做样本\n",
    "\n",
    "t_statistic,p_value=stats.ttest_1samp(a=sample_data,popmean=data.mean() + 5000)\n",
    "print('从样本构造的t统计量 = ',t_statistic)    #t的绝对值越大，p值越小，可参见下面的t-分布图\n",
    "print('单样本双边检验的 p = ',p_value)\n",
    "\n",
    "#画一下样本容量为80时的t分布曲线\n",
    "df=80-1\n",
    "\n",
    "x = np.linspace(stats.t.ppf(0.00000001,df),stats.t.ppf(0.99999999,df),80)  #ppf函数是CDF的逆函数，用来求分位点\n",
    "plt.plot(x, t.pdf(x, df))  #pdf生成概率密度函数表\n",
    "plt.plot((t_statistic,t_statistic),(-0.01,0.4),'-.r')\n",
    "str_legend=('t distribution','calculated t')\n",
    "plt.legend(str_legend)\n",
    "\n",
    "\n",
    "x_025=stats.t.ppf(0.025,df) #注意：同样选取alpha=0.5，做双边test时要给左、右两侧各分配0.25\n",
    "plt.plot((x_025,x_025),(-0.01,0.4),'--g') #这根线是双边检验的左阈值\n",
    "x_975=stats.t.ppf(0.975,df)\n",
    "plt.plot((x_975,x_975),(-0.01,0.4),'--g') #这根线是双边检验的右阈值\n",
    "#print(x_975)\n",
    "\n",
    "rect1=plt.Rectangle((x_025,0),x_975*2,0.5,color='g',alpha=0.25)\n",
    "plt.gca().add_patch(rect1)\n",
    "\n",
    "rect2=plt.Rectangle((x_975,0),5,0.5,color='r',alpha=0.25)\n",
    "plt.gca().add_patch(rect2)\n",
    "rect3=plt.Rectangle((x_025,0),-5,0.5,color='r',alpha=0.25)\n",
    "plt.gca().add_patch(rect3)\n",
    "\n",
    "#plt.ylim(-0.02,0.5)\n",
    "plt.text(x_025,0.4,'accept region')   #图片上不支持中文\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p =  0.00470846938730825 \n",
      "\n",
      "\u001b[1;31m 拒绝 \u001b[0m H0,即总体的均值不等于 13421.026851851851 ，此时错误拒绝H0的概率为 0.00470846938730825 小于显著性水平\n"
     ]
    }
   ],
   "source": [
    "print('p = ',p_value,'\\n')\n",
    "alpha=0.05\n",
    "if p_value>=0.05:\n",
    "    print('\\033[1;31m 接受 \\033[0m H0')\n",
    "else:\n",
    "    print('\\033[1;31m 拒绝 \\033[0m H0,即总体的均值不等于',data.mean() + 5000,'，此时错误拒绝H0的概率为',p_value,'小于显著性水平')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "由显著性检验的原理可知，如果假设H0正确，那么他被错误拒绝的概率不应该大于$\\alpha$的值，即0.5。     \n",
    "而不论是从计算得出的p值还是从自由度79的t-分布的图像中，都可以很容易的看出，实际数据与理论假设H0的偏离达到了“显著”的程度。    \n",
    "由中心极限定理，可以认为几乎不可能出现错误拒绝的情况，因此认为假设是错误的。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 双样本检验\n",
    "利用双样本检验可以将样本集A的统计均值与样本集B的统计均值比较，比如，随机抽取两组样本数据，判断两组样本来源的总体均值是否相等/相近：    \n",
    "从同一个数据集中抽取两次数据，从理论上来说，应该大概率可以接受来源相等的假设："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "第 0 组样本的均值为 7762.225\n",
      "第 1 组样本的均值为 8925.946249999999\n"
     ]
    }
   ],
   "source": [
    "my_sample={}\n",
    "for n in range(2):\n",
    "    my_sample[n]=np.random.choice(a=data,size=80) #从中随机抽取80个做样本\n",
    "    print('第',n,'组样本的均值为',my_sample[n].mean())   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "t =  -0.66244712661189\n",
      "p =  0.5096127404838778\n",
      "\u001b[1;31m 接受 \u001b[0m the H0: 两样本来源的总体均值相等\n"
     ]
    }
   ],
   "source": [
    "#能否下结论说两组样本来源的总体其均值不相等呢？\n",
    "#需进行two sample t-test\n",
    "#H0: 两样本均值相等\n",
    "alpha=0.01   #设置显著性水平\n",
    "t_statistic,p_value=stats.ttest_rel(a=my_sample[0],b=my_sample[1]) \n",
    "#sample size相等的双样本均值比较\n",
    "print('t = ',t_statistic)\n",
    "print('p = ',p_value)\n",
    "if p_value<=alpha:\n",
    "    print('\\033[1;31m 拒绝 \\033[0m H0：两样本来源的总体均值相等 ')\n",
    "else:\n",
    "    print('\\033[1;31m 接受 \\033[0m the H0: 两样本来源的总体均值相等')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "当然，也可以重新抽取不同的年份作为数据集，判断一下是否随机到了相同的年份："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "第1组样本的统计区间为 1978-1982 统计均值为 3934.2737499999994\n",
      "第2组样本的统计区间为 2013-2017 统计均值为 15426.97375\n",
      "t =  -4.901520811433556\n",
      "p =  4.974950696470544e-06\n",
      "\u001b[1;31m 拒绝 \u001b[0m H0(两样本来源的总体均值相等)\n",
      "应该来自不同的年份 \n"
     ]
    }
   ],
   "source": [
    "temp_year_1 = random.randint(0,11)\n",
    "temp_year_2 = random.randint(0,11)\n",
    "data_1 = data_gdp[time_periods[temp_year_1]]\n",
    "data_2 = data_gdp[time_periods[temp_year_2]]\n",
    "new_sample={}\n",
    "my_sample[0]=np.random.choice(a=data_1,size=80) #从中随机抽取80个做样本\n",
    "my_sample[1]=np.random.choice(a=data_2,size=80) #从中随机抽取80个做样本\n",
    "print('第1组样本的统计区间为',time_periods[temp_year_1],'统计均值为',my_sample[0].mean())   \n",
    "print('第2组样本的统计区间为',time_periods[temp_year_2],'统计均值为',my_sample[1].mean())   \n",
    "alpha=0.01   #设置显著性水平\n",
    "t_statistic,p_value=stats.ttest_rel(a=my_sample[0],b=my_sample[1]) \n",
    "#sample size相等的双样本均值比较\n",
    "print('t = ',t_statistic)\n",
    "print('p = ',p_value)\n",
    "if p_value<=alpha:\n",
    "    print('\\033[1;31m 拒绝 \\033[0m H0(两样本来源的总体均值相等)\\n应该来自不同的年份 ')\n",
    "else:\n",
    "    print('\\033[1;31m 接受 \\033[0m the H0(两样本来源的总体均值相等)\\n可能来自相同的年份')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 卡方检验（随堂练习：Titanic-2）\n",
    "数据在各类上的分布是否与给定的概率分布一致？    \n",
    "舱位等级对旅客的存活率有影响吗？仿照按照性别划分的课堂demo，有如下判断："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "各个舱位的频率为 [0.24242424 0.20650954 0.55106622]\n"
     ]
    }
   ],
   "source": [
    "mask1=titanic['Pclass']==1\n",
    "mask2=titanic['Pclass']==2\n",
    "mask3=titanic['Pclass']==3\n",
    "p=np.array([sum(mask1)/(sum(mask1)+sum(mask2)+sum(mask3)),sum(mask2)/(sum(mask1)+sum(mask2)+sum(mask3)),sum(mask3)/(sum(mask1)+sum(mask2)+sum(mask3))])\n",
    "print('各个舱位的频率为',p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "总体存活人数为 342\n"
     ]
    }
   ],
   "source": [
    "mask_survived=titanic['Survived']==1\n",
    "my_population=titanic.loc[mask_survived,'Pclass']\n",
    "# print(type(my_population))\n",
    "pop_size=my_population.count()\n",
    "print('总体存活人数为',pop_size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "如果每个舱位获救机会均等，那么三个舱位分别的存活人数是 [ 82.90909091  70.62626263 188.46464646]\n"
     ]
    }
   ],
   "source": [
    "sample=my_population\n",
    "E=pop_size*p\n",
    "print('如果每个舱位获救机会均等，那么三个舱位分别的存活人数是',E)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "实际的存活人数是  [136  87 119]\n"
     ]
    }
   ],
   "source": [
    "mask1=sample==1\n",
    "mask2=sample==2\n",
    "mask3=sample==3\n",
    "my_set1=sample[mask1]\n",
    "my_set2=sample[mask2]\n",
    "my_set3=sample[mask3]\n",
    "O=np.array([len(my_set1),len(my_set2),len(my_set3)])\n",
    "print('实际的存活人数是 ',O)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "t= 63.39624559772318 \n",
      "p= 1.7126955990497913e-14\n",
      "我们 \u001b[1;31m 拒绝 \u001b[0m 各个舱位具有相同死亡率的假设.\n"
     ]
    }
   ],
   "source": [
    "chi_squard,p_value=stats.chisquare(f_obs=O,f_exp=E)\n",
    "print('t=',chi_squard,'\\np=',p_value)\n",
    "\n",
    "a=0.05\n",
    "df=2\n",
    "if p_value<=a:\n",
    "    print('我们 \\033[1;31m 拒绝 \\033[0m 各个舱位具有相同死亡率的假设.')\n",
    "else:\n",
    "    print('我们 \\033[1;31m 接受 \\033[0m 各个舱位具有相同死亡率的假设')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 卡方检验\n",
    "最后，回到AQUASTAT数据集，我们来判断一下不同的区域的经济平均发展水平是否一致吧（虽然用腿想都知道答案）      \n",
    "首先，输出包含的区域，这里借鉴了[pycon-2017-eda-tutorial](https://github.com/cmawer/pycon-2017-eda-tutorial)的做法，使用一个字典将区域重新划分为简化版的地理区域："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['World | Asia',\n",
       "       'Americas | Central America and Caribbean | Central America',\n",
       "       'Americas | Central America and Caribbean | Greater Antilles',\n",
       "       'Americas | Central America and Caribbean | Lesser Antilles and Bahamas',\n",
       "       'Americas | Northern America | Northern America',\n",
       "       'Americas | Northern America | Mexico',\n",
       "       'Americas | Southern America | Guyana',\n",
       "       'Americas | Southern America | Andean',\n",
       "       'Americas | Southern America | Brazil',\n",
       "       'Americas | Southern America | Southern America', 'World | Africa',\n",
       "       'World | Europe', 'World | Oceania'], dtype=object)"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_data.region.unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "simple_regions ={\n",
    "    'World | Asia':'Asia',\n",
    "    'Americas | Central America and Caribbean | Central America': 'North America',\n",
    "    'Americas | Central America and Caribbean | Greater Antilles': 'North America',\n",
    "    'Americas | Central America and Caribbean | Lesser Antilles and Bahamas': 'North America',\n",
    "    'Americas | Northern America | Northern America': 'North America',\n",
    "    'Americas | Northern America | Mexico': 'North America',\n",
    "    'Americas | Southern America | Guyana':'South America',\n",
    "    'Americas | Southern America | Andean':'South America',\n",
    "    'Americas | Southern America | Brazil':'South America',\n",
    "    'Americas | Southern America | Southern America':'South America', \n",
    "    'World | Africa':'Africa',\n",
    "    'World | Europe':'Europe', \n",
    "    'World | Oceania':'Oceania'\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['Asia' 'North America' 'South America' 'Africa' 'Europe' 'Oceania']\n"
     ]
    }
   ],
   "source": [
    "data_new_region = pd.read_csv('aquastat.csv.gzip', compression='gzip') # 好像会和前面起冲突，就重新读取一遍吧\n",
    "data_new_region.region = data_new_region.region.apply(lambda x: simple_regions[x])\n",
    "new_region = data_new_region.region.unique()\n",
    "print(new_region)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "然后，将不同区域的数据储存到不同的DataFrame中，然后把人均GDP数据切分出来，最后随机出一个特定年份，对该统计区间内的数据进行操作："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "def subregion(data, region):\n",
    "    return data[data.region==region]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_div_region={}\n",
    "data_div_region_gdp={}\n",
    "for n in range(6):\n",
    "    data_div_region[n] = subregion(data_new_region, new_region[n])\n",
    "    data_div_region_gdp[n] = slice_by_variable(data_div_region[n], 'gdp_per_capita')\n",
    "    data_div_region_gdp[n] = data_div_region_gdp[n].dropna(axis=0)\n",
    "# print(data_div_region_gdp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[1;31m统计区间： 1998-2002 \u001b[0m\n"
     ]
    }
   ],
   "source": [
    "data_final={}\n",
    "temp_year = random.randint(0,11)\n",
    "print('\\033[1;31m统计区间：',time_periods[temp_year],'\\033[0m')\n",
    "for n in range(6):\n",
    "    data_final[n] = data_div_region_gdp[n][time_periods[temp_year]]\n",
    "# print(data_final)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "由于统计数据和统计方便的缘故，我们设定一个“发达指标”如下：     \n",
    "<u>如果一个国家的人均GDP超过了全球所有国家的中位数，那我们就认为这个国家是“发达国家”，如果全球发展水平一致的话，那么所有地区的“发达国家”占比都应该在50%左右。</u>  \n",
    "同时，由于大洋洲就三个国家，容易引起极端情况，因此就不再考虑大洋洲国家的数据。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[1;31m对于统计区间： 1998-2002 ,其全球人均GDP中位数为： 2079.0 \u001b[0m\n"
     ]
    }
   ],
   "source": [
    "mid = slice_by_variable(my_data, 'gdp_per_capita')[time_periods[temp_year]].median()\n",
    "print('\\033[1;31m对于统计区间：',time_periods[temp_year],',其全球人均GDP中位数为：',mid,'\\033[0m')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "首先统计不同地区的国家占总国家数量的频率："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "各个地区的频率为 [0.19047619 0.18095238 0.1047619  0.35238095 0.17142857]\n"
     ]
    }
   ],
   "source": [
    "sum_all = len(data_final[0])+len(data_final[1])+len(data_final[2])+len(data_final[3])+len(data_final[4])\n",
    "p=np.array([len(data_final[0])/sum_all,len(data_final[1])/sum_all,len(data_final[2])/sum_all,len(data_final[3])/sum_all,len(data_final[4])/sum_all])\n",
    "print('各个地区的频率为',p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "然后，根据“发达国家”的定义，统计出各个地区的发达国家的数量与全球所有发达国家的数量："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "for n in range(5):\n",
    "    data_final[n] = data_final[n][data_final[n]>=mid]\n",
    "# print(data_final)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "sum_all = len(data_final[0])+len(data_final[1])+len(data_final[2])+len(data_final[3])+len(data_final[4])\n",
    "q=np.array([len(data_final[0])/sum_all,len(data_final[1])/sum_all,len(data_final[2])/sum_all,len(data_final[3])/sum_all,len(data_final[4])/sum_all])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`E`代表猜想全球发展水平一致的情况下，各个地区应该有的发达国家的数量      \n",
    "`O`代表实际情况下，各个地区的发达国家的数量"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[11.23809524 10.67619048  6.18095238 20.79047619 10.11428571]\n",
      "[10. 15.  8.  8. 18.]\n"
     ]
    }
   ],
   "source": [
    "E = p * sum_all\n",
    "O = q * sum_all\n",
    "print(E)\n",
    "print(O)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "代入检验公式，检验假设是否成立"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "t= 16.43986042647952 \n",
      "p= 0.0024823173376233194\n",
      "我们 \u001b[1;31m 拒绝 \u001b[0m 各个地区“发达程度”相同的假设\n"
     ]
    }
   ],
   "source": [
    "chi_squard,p_value=stats.chisquare(f_obs=O,f_exp=E)\n",
    "print('t=',chi_squard,'\\np=',p_value)\n",
    "\n",
    "a=0.05\n",
    "df=4\n",
    "if p_value<=a:\n",
    "    print('我们 \\033[1;31m 拒绝 \\033[0m 各个地区“发达程度”相同的假设')\n",
    "else:\n",
    "    print('我们 \\033[1;31m 接受 \\033[0m 各个地区“发达程度”相同的假设')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "当然，这种判断方法是很粗糙的，有很多没有考虑到的因素，比如说，是否应该根据人口数量进行加权，计算出各个地区根据人口分布应该有的“发达国家”数量，而非简单地按照国家数量计算。    \n",
    "同时，简单地按照人均GDP的中位数进行判断是否准确也是一个问题。     \n",
    "（不过结论倒是对的，毕竟不用判断也知道这个问题的答案）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 总结\n",
    "通过这次作业学习了区间估计和假设检验的方法，练习了如何通过数据科学的方法验证自己对大数据的猜想。     \n",
    "卡方检验也可以完成对数据关联性的一些简单检验，比如使用各个国家的水资源数据进行加权运算，可以得到“发达国家数”按照水资源多寡的分布，再与实际的发达国家进行对比，就可以验证国家“发达”与水资源多寡之间的关联。"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
